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I. INTRODUCTION 


An increasingly important measure of the usefulness of aircraft flight test 
data is its ability to estimate stability and control derivatives . Improvements 
in the algorithms to accomplish this parameter identification task have advanced 
to the point where the choice of control inputs may be a limiting factor in the 
attainable accuracy of these estimated stability and control derivatives. An input 
design requirement, therefore, has arisen because of the need to improve the 
efficiency of flight testing by obtaining more accurate estimates from response data 
in less time , 

The importance of choosing appropriate control inputs and exciting specific 
aircraft modes for extracting stability and control derivatives from aircraft flight 
testing has long been recognized. As early as 1951, Milliken^ summarized the 
studies conducted in defining good input signals in this statement: 11 It would appear 

that an optimal input in a given case is that which best excites the frequency of 
interest, and, hence, its (the input signal) harmonic content should be examined 
before the test to ensure that it is suitable" . Good inputs could resolve parameter 
identifiability problems and improve confidences on estimates of stability and control 
derivatives obtained from the resulting flight test data. In other words, with spe- 
cially chosen inputs, the same accuracy on parameter estimates can be obtained in 
much shorter flight test time than with conventional inputs . Shorter flight tests can 
lead to a saving in time required for stability and control testing and the computation 
requirements for extraction of aerodynamic derivatives. In addition, these inputs 
can be chosen specifically to satisfy the ultimate flight test objective such as control 
systems design, simulator parameter specifications, response prediction, aero- 
dynamic model validation, or handling qualities evaluation. 

There are many factors that must be considered when choosing inputs for 
flight tests . These include: 

(a) Pilot Acceptability . If the flight test is to be carried out with a pilot 

onboard the aircraft, it is necessary that the control inputs be accept- 
able to the pilot. The inputs should not maneuver the aircraft into a 
flight region from which a pilot cannot recover. In addition, the inputs 


should be reproducible by the pilot . 


(b) Instrumentation . The inputs must consider specific instruments avail- 
able, and their dynamic range and accuracy. The primary impact of 
the instruments on input design is on the signal/noise ratios which the 
response must have for sufficiently accurate data. 

(c) Parameter Identification Technique . Many parameter extraction 
methods require a certain class of inputs (e.g. , sinusoidal inputs for 
transfer function identification, random inputs for correlation tech- 
niques) . More advanced techniques of parameter identification tend 
not to rely on such specific classes of inputs, but do require inputs 
which maximize some function of the sensitivity of the output responses 
to parameters . 

(d) Modeling Assumptions . The inputs that are designed must also consid- 
er the model that is assumed. For example, inputs chosen for a linear 
mode should not cause such large aircraft motions that the assumption 
of constant stability and control derivatives is invalid . 

(e) Aircraft Structural Constraints . The aircraft maneuvers produced by 
the inputs should not cause the structural loads to increase beyond the 
design stresses of various aircraft components . 

(f) Objective of Parameter Identification . This is one of the most important 
and least understood of input design requirements . It is now known - 
that there may be a significant difference in inputs which allow more 
accurate estimates of parameters for control system design as opposed 
to those inputs required for estimates of handling quality coefficients . 
Unfortunately, there Eire not extensive systematic techniques for relat- 
ing the input design to the identification objective . 

(g) Output Sensitivity . Measured aircraft response resulting from these 
inputs should be most sensitive to the parameters of interest and less 
sensitive to other, possibly unknown, parameters of interest. 
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Figure 1. 1 illustrates the use of specially designed inputs for flight test de- 
sign. This procedure is based on specification of overall aircraft characteristics, 
instrumentation, and parameter identification objectives . The use of an algorithm 
to design the input may preclude consideration of other possible constraints (due to 
computational complexity) so some iteration to meet other constraints not considered 
may be desirable. Once flight tests are completed, and parameter extraction per- 
formed, the identification results may be used to design other inputs to further 
improve accuracy or as a priori estimates to be used in identifying other parameters 
not originally considered . 

[2 3 ] 

The present work is an extension of a previous study in which a basic 

approach to input design was formulated based on optimization of a function of the 
sensitivity of the aircraft response to the aerodynamic derivations (e.g . , the infor- 
mation matrix, M) . Though this initial effort did establish the feasibility of the 
approach, its application to flight test requirements was difficult. To facilitate this 
application, the present effort was initiated. 

The basic objective was to extend the formulation of input .design procedure 
to make it more useful for flight test application and to meet the requirements listed 
above. This task was to be achieved by extension of the original time domain 
method and also by development of a frequency domain technique. In the course of 
the work, it was realized that the original approach could be reformulated by opti- 
mizing a function of the inverse of the information matrix (i .e . , the dispersion 
matrix, D) for a significant expansion of capability for the resulting inputs to meet 
flight test requirements . Implementation of this reformulation in both the time and 
frequency domain demonstrated both the computational feasibility and the desired 
improvement. Subsequent computation revealed that the time domain method is 
more applicable to most flight test objectives and that the frequency domain approach 
was more useful for test conditions where a steady state condition would be estab- 
lished. 

The capabilities afforded by use of these two techniques includes the follow- 
ing: 
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Figure 1. 1 Use of Optimal Inputs for Flight Test Design 














(a) 


Choice of Criteria In Choosing Inputs: Given that the constraints 

listed above are satisfied, there is still flexibility in selecting certain 
input objectives . Specifically, it may be necessary that some para- 
meters must be very accurately identified, while others may not be so.* 
Such an objective is best met by using a weighted trace of the informa- 
tion or dispersion matrix. Alternately, the eventual use of the estimates 
may not be clearly known, but it is desired to maximize the overall 
identifiability of the responses for whatever purpose (e.g., handling 
quality evaluation, control system design, etc.). Such an objective 
may be achieved by using the determinant of the information or disper- 
sion matrix. The two techniques of this work allow use of any of these 
options (e .g . , Tr M, Tr D , |Mj , |D| ) . It will be shown that choosing 
inputs based on the trace of the dispersion matrix option should give 
the most accurate estimates, in general. 

(b) Ability To Be Able To Identify a Large Number of Parameters : Pre- 

viously, input design techniques for linear systems were limited by 
computational requirements to low order aircraft models involving only 
a few parameters. The techniques developed for this work design 
inputs for a significantly larger number of parameters to be identified. 
This capability allows compression of flight test time to acquire the 
data to estimate the most derivatives . 

These capabilities have been evaluated on simulated aircraft data of the 
Buffalo C-8 longitudinal response and the Lockheed Jet Star lateral response . The 
methods of evaluation are: 

(1) Comparison of the two input design approaches between each other 
and also against pulse and doublet-type inputs . 

* __ 

Alternately, some derivatives, such as , may be well known from previous 

wind tunnel or flight tests, and others, such as Cg , may not be confidently estab- 

P 

lished. Inputs could then be designed only for the latter group. This situation is 
not emphasized in this work. 
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(2) Comparison of the input designs against such factors as ease of imple- 
mentation and levels of aircraft response . 

The method of comparison is based on the standard deviations of the estimates from 
one input versus that of another . These criterion values are obtained as the square 
roots of the diagonal of the dispersion matrix corresponding to the input under 
evaluation . 

The organization of this report is as follows. Chapter II presents a review of 
the developments of input design for parameter extraction. Chapter III, together 
with Appendices A, B, C, and D, discuss in detail the theoretical background of the 
time domain and the frequency domain methods . Chapter IV presents numerical 
results on optimal flight control inputs obtained using the time domain technique . 
Similar results, for longitudinal motions of a C-8 aircraft and lateral motions of a 
Jet Star, are given in Chapter V . Chapter VI evaluates the inputs under off-design 
conditions and by approximating the inputs by a series of steps . The results and 
conclusions are summarized in Chapter VII. 
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n. REVIEW OF INPUT DESIGN TECHNOLOGY 

2 . 1 CONVENTIONAL INPUTS IN AIRCRAFT APPLICATIONS 

Since the first efforts of applying parameter extraction technology to aircraft 

flight test data, many different control inputs have been used. Many flight 

tests are presently aimed at determining natural frequency, damping ratio, etc. , 

[5] 

of a specific mode and steady state gains . Most of the inputs are selected on the 
needs of simple parameter extraction procedures . One commonly used input is the 
frequency sweep . ^ j n this approach, the aircraft is. excited by sinusoidal inputs 
over a range of frequencies , usually around the natural frequency of the mode , 
until a steady state is reached at each frequency. The parameters of a suitable lin- 
ear model are selected to obtain the best fit to the variation with input frequency 
of the output/input amplitude ratio and phase difference . These inputs work 

satisfactorily but require much flight test time. With the development of more 

r 70 1 

sophisticated parameter extraction methods, other inputs have been tried. ’ 

Pulse inputs are used sometimes, and the frequency response is obtained by 
taking the Fourier transform of the output and the input at discrete points . 

These inputs are limited to simple low order linear systems . 

Doublets, steps, and finite duration pulse inputs are generally used to identi- 
fy aircraft parameters in both linear and nonlinear flight regimes. However, the 
estimates of certain parameters may be quite poor and, in some cases, a set of para- 
meters may not be identifiable at all. A possible result is that excessive flight test 
time may be required to get good estimates of all the parameters. Optimal inputs 
consider the identifiability of each stability and control derivative directly. They 
can be tuned to obtain better identifiability of the overall parameter set or tuned to 
identifying only particular parameters of primary interest . 

2.2 OPTIMAL INPUT DESIGN METHODS 

The use of analytical techniques in input design may be considered to have 
been initiated by Fisher 1 who gave a quantitative meaning to the knowledge 
about a certain set of parameters through definition of the information matrix, M. 
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sophisticated parameter extraction methods, other inputs have been tried. * 

Pulse inputs are used sometimes, and the frequency response is obtained by 
taking the Fourier transform of the output and the input at discrete points. 

These inputs are limited to simple low order linear systems. 

Doublets , steps, and finite duration pulse inputs are generally used to identi- 
fy aircraft parameters in both linear and nonlinear flight regimes . However, the 
estimates of certain parameters may be quite poor and, in some cases, a set of para- 
meters may not be identifiable at all. A possible result is that excessive flight test 
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consider the identifiability of each stability and control derivative directly. They 
can be tuned to obtain better identifiability of the overall parameter set or tuned to 
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[Q ] 

been initiated by Fisher 1 who gave a quantitative meaning to the knowledge 
about a certain set of parameters through definition of the information matrix, M. 
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Later, Cramer and Rao^^ showed that the inverse of the information matrix, referred 
to as the dispersion matrix, is abound on parameter error covariances, i.e., 


cov { (0 - 0) (0 - 0) T } > M 1 A D 


( 2 . 1 ) 


where 0 is the actual value of a parameter and 0 is the estimate. This is called 
the Cramer-Rao lower bound. The dispersion matrix is extremely useful because 
it relates to physically meaningful quantities (parameter estimation errors) and 
gives a method of comparing different experiments . Most of the analytical methods 
in input design use a function of information or dispersion matrix as the extremiz- 
ing criterion. It is assumed that an efficient parameter identification algorithm, 
which can extract all available information about the parameters from data, is 
used and the Cramer-Rao lower bound is met with equality. This is important 
since this makes it possible to design inputs independent of the parameter 
extraction procedure . 

2.2.1 Inputs for Regression Systems 

Statistical and optimization concepts for selection of inputs were first used 
for regression experiments (e.g. , static systems) . Suppose it is possible to take 
noisy measurement of linear combinations of a set of unknown parameters,©, i.e. , 

y = F (u)0 + n (2.2) 

where y is the observations and n is random noise. The set of possible linear 

combinations is defined by the control u. The maximum number of observations 

is N . The input design problem consists of finding a set of u's and the number 

of observations at each u to get the "best" estimate of parameters, 0, under the 

above constraints . Some of the earliest work in this field was done by Kiefer and 

Wolfowitz, where they proved the equivalence of two extremum problems. 

They showed that optimizing a certain criterion in parameter space (e.g., sum of 

covariance of parameter estimates) is equivalent to optimizing a certain other 

criterion in sample or output space. This is an important conclusion and resulted 

r 12 1 

in a number of significant contributions, e.g. , Kiefer, and Karlin and Studden. 
An excellent summary of these methods is given in Fedorov. 
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br> 2.2.2 Inputs for Dynamic Systems 


Earliest work in the field of input design for dynamic systems was done by 
Levin and Litman and Huggins. These authors designed inputs for un- 

known parameters in system impulse response based on least square estimation. 

In addition, Litman and Huggins considered an infinite observation time. 

. [17] 

Levadi was able to put the problem in a more general framework. In his 
approach, the observation time is finite and there are constraints on total input 
energy. By using a variational procedure, the trace of the error covariance 
matrix is minimized leading to a nonlinear Fredholm equation. 

[18] 

Aoki and Staley designed inputs for single-input, single-output systems 

based on maximizing the trace of the information matrix. They considered discrete 

time representation of a dynamic system and showed that the energy bounded 

optimal input is the eigenvector of a certain matrix corresponding to its maximum 

[ 19] 

eigenvalue . Nahi and Wallis have also considered this problem but did not 

[ 2 ] 

come up with a general algorithm. Mehra proposed an algorithm which maxi- 
mizes the trace of the information matrix for multi- input, multi-output systems. - 
Several methods for solving the resulting two point boundary value problem were 
also given. One method required solving a Riccati equation and was tried with 
partial success. The designed inputs, using this method, are unsuitable because, 
in general, maximizing the trace of the information matrix does not ensure that 
the covariances on parameter estimates are small. In extreme circumstances, it 
may give a singular information matrix, which gives infinite covariance in certain 
directions in the parameter space. Reid^^ and Goodwin et al. design 
inputs based on the trace of the dispersion matrix (sum of parameter error covari- 
ances). They use direct gradient procedures, which require excessive computa- 
tion time even for simple systems. The method, they propose, is unusable for 
multi- input, multi-output systems and for systems with more than a few unknown 
parameters . 

[ 22 ] 

Mehra , in a novel approach, uses the steady state assumption to convert 
a linear constant coefficient system into its frequency domain representation. He 
demonstrates the procedure for determining optimal input spectra to minimize the 
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determinant, trace or any of a variety of functions of the dispersion matrix. 

[ 23 ] 

Viort also considers a similar problem. 

Until now, inputs have been computed for only a few simple practical sys- 
tems . There have been two major barriers limiting the determination of optimal 
inputs . Most techniques, developed to date, can handle only the trace of the 
information matrix as the optimizing criterion and in many cases these inputs pro- 
duce either marginal or no improvement in parameter estimation accuracy. 
Secondly, the computation time required is so large for high order systems (e.g., 
more than three states) with many unknown parameters (e.g, ten) that it makes 

the actual determination of inputs infeasible. This is because largely brute force 

[2 20 21 ] 

methods have been used in the past. 9 9 In aircraft applications, Stepner 

[31 

and Mehra computed inputs which maximize the trace of the information matrix 

for identifying five parameters in the longitudinal short period mode using the two 

state approximation. The computed input gave better estimates of three parameters 

but poorer estimates of two parameters as compared to a conventional doublet 

[ 24 ] 

input. Swanson and Bellville have used some of these techniques to design 
inputs to identify parameters in certain biological systems . 

In summary, previous input design techniques for dynamic systems have 
demonstrated the potential of improving the capability to identify aerodynamic 
derivatives from flight data. These techniques have, however, been limited in the 
flexibility they allow to meet important flight test requirements such as obtaining 
high accuracy for specific derivatives. In addition, such previous approaches 
have not demonstrated the capability to provide inputs for estimation of a large 
number of derivatives within reasonable computation limits. In the following 
chapter, an input design method is discussed which is directed toward alleviating 
these problems . 
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El . METHOD OF TIME AND FREQUENCY INPUT DESIGN 
FOR DYNAMIC SYSTEMS 


3 . 1 INTRODUCTION 

This chapter describes a method consisting of two different techniques for 
design of input signals which provide estimates of unknown parameters in linear 
time -invariant systems. The first technique uses the time domain representation 
of system dynamics and develops methods to compute the time history of control 
input sequence for any duration of the experiment. In the second technique, the 
system is assumed to be in oscillatory steady state and a frequency domain repre- 
sentation is utilized. This gives the optimal control input spectrum. The corres- 
ponding time history is "optimal" only for long experiments. The time domain 
approach is computationally much more complicated than the frequency domain 
approach. The two approaches are, therefore, complementary . The frequency 
domain approach is suitable for long experiments and the time domain approach 
should be used for short and medium experiments . 

A computation algorithm based on eigenvalue-eigenvector decomposition is 
developed to solve the time domain problem. A new sensitivity functions reduction 
method affords considerable savings in computation time by decreasing the order 
of the problem. These two algorithms have made implementation of the time domain 
algorithm more feasible for practical systems because they allow design of inputs 
for much higher order systems than previously reported. 

3 . 2 PROBLEM STATEMENT 

Consider a linear, time-invariant, dynamic system following the differential 
equation: 


x = Fx + Gu 


x(0) = 0 


where 


0 < t< T 


(3.1) 
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x is an nxl state vector, 
u is a qxl input vector, and 

F and G are appropriate matrices which depend on m unknown 
parameters 9. 

Let there be continuous noisy measurements of p linear combinations of state 
variables . 

y = Hx + v (3.2) 

where 

y is a pxl measurement vector, 

v is a pxl white noise vector with zero mean and power spectral 
density R , and 

H is a pxn matrix which is a function of parameters 0. 

The problem is to choose u from a class of inputs to obtain "best" estimates of 
the unknown parameters . A total energy constraint is imposed on the input to 
limit state and control input excursions. This method of state deviation constraint 
is indirect. The more desirable and direct method of including quadratic penalty 
on state in the cost function is a difficult analytical and computational extension 
to the present approach of limiting total input energy alone, 

r T 

I u^u dt = E . (3.3) 

q • » 


3 . 3 CRITERIA OF OPTIMALITY 

It is usually not possible to find an input which gives better estimates of all 
parameters in a given system than any other input. The optimal input is deter- 
mined by giving suitable importance to different parameters. Let M be the infor- 
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mation matrix for parameters 0 resulting from an input u. Then the dispersion 
matrix is defined as 

DAM -1 (3.4) 

and the Cramer-Rao lower bound (Equation (2.1)) gives 


cov (0) > D (3.5) 

where 0 is the estimated value of 0. In general, it is necessary to make a trade- 
off among the accuracies on estimates of unknown parameters in the system. 

Based on Equation (3.5) , several optimality criteria have been proposed and 
used. Some of these criteria are not altogether appropriate in that they do not 
ensure small estimation errors , but have been used because it is easier to compute 
the corresponding "optimal 11 input. There are three classes of criteria which 
have received special attention in the past. 

(a) Linear functional of the information matrix. 

= max ^(M) (3.6) 

u 

(b) The determinant of the dispersion matrix. 

= min |D| (3 . 7) 

u 

(c) Linear functional of the dispersion matrix. 

J 3 = min «2?(M) (3.8) 

u 

is such that for two positive semi-definite matrices A and B and a constant c 
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(a) ^(A) > 0 

(b) i^(A + B) = £?(A) + i?(B) (3.9) 

(c) i?(cA) = c^(A) 

Examples of linear operator ^ are the trace and the weighted trace. 

maximizes the total or partial sum of information of all the parameters 
or of a linear transformation of parameters . If the linear operator is the trace, 
the total information of all parameters is maximized. This may, however, lead to 
an almost singular information matrix with large terms on the diagonal. Then the 
dispersion matrix, which is a lower bound on parameter error covariances, has 
large diagonal terms. Therefore, this optimality criterion is not very suitable. 

It is used mainly because of its simplicity. 

The dispersion matrix which is positive-definite, in the light of the Cramer- 
Rao lower bound, can be looked upon as a hyper ellipsoid of uncertainty in the 
parameter space . works with the determinant of the dispersion matrix and 
minimizes the volume of the uncertainty ellipsoid. 

minimizes a weighted sum of covariances of parameter estimates (or some 
linear combinations of parameters) . The weighting matrix serves two purposes . 
Since the covariances of different parameters have different units, it converts each 
term in the sum to the same units . Secondly, the weighting matrix offers a tremen- 
dous flexibility because it is possible to assign varying importance to parameters, 
through weights on their nondimensional covariance. This is considered to be 
one of the most suitable performance criteria since it works with parameter esti- 
mate covariances directly. 

In the next sections, we indicate how these different criteria can be handled 
in the time domain approach and the frequency domain approach. 
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3.4 TIME DOMAIN INPUT DESIGN TECHNIQUE FOR DYNAMIC SYSTEMS 


In the past it has been possible to work only with a linear function of the 
information matrix, in particular, the trace. Under the present effort, methods 
have been developed which make it possible to optimize the determinant or weighted 
trace of the dispersion matrix. The details of the theory behind these methods is 
given in Appendix A. Here, we describe the algorithms and then indicate numeri- 
cal procedures which are used to solve the resulting equations. 


yhted Trace of the Information Matrix 


It is shown in Appendix A that maximizing the trace of the information leads 

to an eigenvalue problem of a positive self adjoint function. It is possible to re- 

[ 2 ] 

formulate it as a two-point boundary value problem. 



L H e R e lH e 



x 0 (O) = 0 X(T) = 0 

u op t = - ^ G e *• 


C3.10) 

(3.11) 


The smallest value of constant [i for a nontrivial solution to the two-point boundary 
value problem gives the optimal control input. The matrices Xq, Fq, G^, H^ and 
Rq acre 
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( 3 . 12 ) 


( 3 . 13 ) 


( 3 . 14 ) 


( 3 . 15 ) 
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This two point boundary value problem is solved as described in Section 
3.4.3 and in Appendix C . 

A weighted trace of the information matrix is maximized by defining a new 
set of m parameters <p related to 0 as 


0 = Cep 


(3.16) 


The information matrices for cp and 0 are related to each other as 

M <p = cX M e c 


(3.17) 


Therefore , 


Tr(M ) =Tr(C T M 0 C) 
= Tr(CC T M 0 ) 


(3.18) 


Thus, maximizing the trace of Mcp is equivalent to maximizing a weighted trace 
of the information matrix for parameters 0. The sensitivities of the state vector to 
parameters tp and 0 are related by the following transformation: 



(3.19) 


It is clear that the parameter transformation of Equation (3.16) enables us to maximize 
a weighted trace of the information matrix, when the weighting matrix is symmetric 
and positive semi-definite . The symmetry is not a restriction and positive semi- 
definiteness is required in the light of condition (a). Equation (3.9). 


3.4.2 Determinant or Weighted Trace of the Dispersion Matrix 

The idea behind minimizing the determinant or weighted trace of the disper- 
sion matrix is presented in Appendix A, Sections A . 4 and A . 5 . We present here 
the algorithm used in the computation of optimal input. It is an iterative pro- 
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cedure with convergence to a stationary point. 


Algorithm (see flowchart in Figure 3.1): 


(1) Choose any input u(t) with energy E which gives a nonsingular in- 
formation matrix,M o . 

(2) Find an input u m (t) with energy E to maximize cp(u), such that 


r 

r. 


u (t) u (t) dt > 0 
o m — 


where , 


cp = Tr(M ^M) to minimize |D| 


= Tr(WM Q ^M M q ^) to minimize Tr(WD) 


(3.20) 

(3.21) 


Both these criteria can be recast as maximizing a weighted trace of 
the information matrix. The matrix C of Equations (3. 16) to (3.19) is 

C = M to minimize | D| (3.22) 


1/2 - 1 

= W ' M q to minimize Tr(WD) 

(3) The information matrix for input au Q (t) + Pu m (t) is 


(3.23) 


M, = a 2 M + B 2 M + 2aBM (3.24) 

1 o r m r om 

where M is the information matrix for input u (t) and M is a 
m m om 

"cross information matrix" for inputs u (t) and u (t) and is defined 

o m 

in Appendix A . The energy constraint on the input requires 
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(3.25) 


a 2 + p 2 + 2ap 

o 

Use Equations (3.24) and (3.25) to find a p between 0 and 1 which 
optimizes the criterion function for M.^. It is shown in Appendix A that 
if the input u Q (t) is not optimal, it is always possible to bring about an 
improvement in the performance index for this choice of p. 

(4) Check to see if the termination criterion is met. One of the following 
can be used . 

(a) The information matrix does not change substantially from one 
step to the next, or if the optimizing function is not improving 
significantly . 

(b) The value of P which optimizes the value of the desired func- 
tion is approaching zero. In other words, very little power is 
being placed at newly chosen frequencies . 

(c) The maximum value of the function $ is not much higher than 
the maximum value for the optimum design (i.e. , m to minimize 
|D| and Tr(VJD ) to minimize TrCWD)) . If the termination 
criterion is not met, repeat from step (2) . 

(5) Check if the design is globally optimum. 

It is clear that, in this technique it is necessary to maximize a 
weighted trace of the information matrix in each iteration. 

Unless otherwise mentioned, in all computations reported in Chapter IV, 
trace or weighted trace of the dispersion matrix is used as the optimality 
criterion . This has proved to be a more useful criterion in general than the trace 
or weighted trace of the information matrix. 
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3.4.3 Sensitivity Functions Reduction 


The two point boundary value problem requires a solution of n(m+l) sensi- 
tivity equations and n(m+l) adjoint equations . It is shown in Appendix D that the 
state sensitivities to system parameters are not independent of each other, in 
general. The number of sensitivity equations can, therefore, be reduced by a 
proper linear mapping and propagating only the independent equations . The 
maximum number of equations to obtain all sensitivity functions is n(q+l) . In 
many practical cases , it is smaller . 

The theory of how the sensitivity functions reductions can be carried out 
[25] 

is given by Gupta and in Appendix D . It is necessary to work only with the 
controllable part of (Fq, Gq) to obtain all sensitivities (Fq and Gq are defined in 
(3.15)). The uncontrollable subspace of (Fq, Gq) is dropped. Also, the states 
unobservable through (Fq, Hq) do not affect the performance index (i.e., the 
trace of the information matrix) . For example, when all parameters are in G, 
the system states, x, are unobservable through Hq and can be discarded. This 
happens rarely and is not incorporated in the algorithm. 

The implemented algorithm is given here and the flowchart is illustrated in 
Figure 3.2. 

1 . 

2 . 


3. 

2k. 1 < n + k. 
i i 


The linearly dependent columns in G are merged. Then the structur- 
ally uncontrollable states in (F, G) are dropped. 

Matrices Fq and Gq are formed. k^,k£, . . . ,k^ of Equation (D.22) are 
determined and are used to choose (q+l)n appropriate columns from 
the controllability matrix of (Fq,Gq) . 


The dimension k!^ of state space controllable from each input u^ alone 
is determined. If for any i 
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only the first 2kj columns involving in the controllability matrix 
are considered, 

4. The remaining columns are checked for linear independence. The 
Gram-Schmidt procedure is used to drop columns, which are linearly 
dependent on other columns . The set of remaining columns is Q-^. 

5. Any pseudo-inverse of is determined. Equation (D . 7) is used to 

compute F and G . 
r c c 

6. Equation (D.6) is solved for x^(t) and Equation (D.8) is used to find 
x 0 (t) at the desired points . 

3.4.4 Solution to the Two Point Boundary Value Problem 


Several solution techniques have been suggested to solve the two point 
boundary value problem of Equation (3. 10) . The Riccati equation method suggest- 

f 2 1 

ed by Mehra has been tried with limited success . In this method, a p is chosen 
and a certain Riccati equation solved to determine the experiment duration where 
the elements of the Riccati matrix become large (theoretically infinite) . The para- 
meter p corresponding to the desired experiment time is determined iteratively . 
The problem with this method is that it is difficult to determine numerically the 
time at which the elements become large. Also, it usually does not give good in- 
sight into the nature of optimal inputs . 


A new method, which uses the symplectic properties of the Hamiltonian 
matrix <7? of the two point boundary value problem, has been developed and is 
described in Appendix C. The eigenvalues of the Hamiltonian occur in pairs 
3P+ and Let the corresponding eigenvector matrix be 



(3,26) 
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with normalization 


T T 

A_ X + - X_ A + = I 


(3.27) 


Then the two point boundary value problem has a nontrivial solution if 


u = a; 1 a. 


-sp.t Y -i Y 

e + X X , 




has at least one eigenvalue equal to one. 


(3.28) 


A computer program has been written to solve Equation (3.10) using eigen- 
vector decomposition. It consists of the following steps: 


(1) A reasonable U is chosen. It can be shown that 

o 


0J 1 
9E | x 


(3.29) 


And since p does not depend upon energy E in u , 

Tr(M)=E/n ! (3.30) 

So can be selected from a knowledge of the energy and the expect- 
ed value of Tr(M) . Alternatively , one could find the value of the 
performance index for a reasonable input which satisfies the energy 
constraint and then apply a suitable correction factor to | x determined 
using Equation (3.29). Choosing a good initial |a Q is important in 
obtaining quick convergence . 


(2) Eigenvalues and eigenvectors of the Hamiltonian are determined and 
A + ^A_ and X_ X + are computed. Starting from T = 0, eigenv alues of 
U are determined for increasing T in steps of AT until a point T q 
is reached where one eigenvalue of U is "close" to 1. 
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ure 3.3 Solution of Two Point Boundary Value Problem Using 
Eigenvalue-Eigenvector Decomposition 
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(3) p is updated for small changes, AT, until the desired T is reached. 
Thereafter, a correction in 1 1 is applied to bring the eigenvalue 
very close to one . 

(4) » the r ight eigenvector of U corresponding to the eigenvalue 
close to one, and X(0) are determined. The states Xq( 0 and input 
u(t) are obtained from the expressions 

x 0 (t) = -X + y R ^ +X_ e~^+ t xjx(0) (3.31) 

u(t) = n(G 0 A + e'^ (T ~ t) y£ 1) - gJ A_ e~ *+* xfx(O)} (3.32) 

(5) The input energy and information matrix are computed. Since the 
system is linear, the states, inputs, and information matrix can be 
scaled for desired energy in control inputs . 

Remarks 

1. In most cases which have been tried, only a few eigenvalues of the 
Hamiltonian are oscillatory or have low damping. Other eigenvalues 
can usually be discarded since they give additional eigenvalues of 

U close to zero. This reduces the size of U. If all but one eigen- 
value of are highly damped, it is usually possible to compute T 

for a \i through knowledge of and S, by hand calculation. 

2. X(0) can also be computed using the following expression: 

U0)=e-^ T X-. T y 1R (3.33) 

It is easier to use this expression since the matrices on the right hand 
side are available in the program. 

3. Equation (3.32) can be simplified to give analytic expressions for 
optimal inputs . 
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3.5 FREQUENCY DOMAIN INPUT DESIGN 


If the duration of the experiment is much longer than the system character- 
istic time, it is possible to design inputs based on a variety of criteria quickly by 

[22] 

making the assumption of steady state. In an ingenious approach, Mehra 
converts a linear time-invariant system into its frequency domain representation. 
This eliminates the dynamics of the system. The parameter estimation problem 
becomes a regression problem in which input frequencies and the power in each 
frequency are the control parameters , These parameters, which define the input 
design, are chosen by an iterative procedure. 

Consider the state space representation of a discrete-time linear system 

x(k + 1) = cpx(k) + Gu(k) k=l, 2, 3, . . . ,N (3.34) 


and the noisy measurements 


y(k) = Hx(k) + v(k) 


(3.35) 


(p , G,and H are appropriate matrices and contain m unknown parameters 0, 
Fourier transform (3.34) and (3.35) to get 


27tn 

y(n) = H(e ^ N - (p) G u(n) + v(n) 
A T (n , 0) u(n) + v(n) 


(3.36) 


As the number of sample points increases, the information matrix per sample 
approaches 


M = 


r 71 * 

1 O / 9T -1, . 3T 

27t Re J ae s w (to) ae 


dF (oo) 
uu 


-7T 


(3.37) 
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where F is the spectral distribution function of u and S is the spectral 
uu vv r 

density of v and superscript '*■ denotes the conjugate transpose of a matrix. 


Based on Appendix B , the following algorithm can be shown to converge to 
the optimum input. See Figure 3.4 for flowchart. 


(1) Choose a nondegenerate input f (co) (i.e., consisting of more than 


frequencies , with a finite power in each frequency) 


(2) Compute the function v|/ (co , f ) and find the value of to where it is 

maximum . Call it co , where 

o 

4-«o.f o ) =Re [ttS^C®) $5^ D« o ) (3.38) 

to minimize |D| 

and 

*<“.0 = 2 [d« 0 ) to) D(f o )j (3.39) 

to minimize a linear function SB of D . 


(3) Evaluate the normalized information matrix at co . Call it M(co ) . 

o o 

(4) Update the design 

f = (1 - a ) f + a f(co ) 0 < a < 1 (3.40) 

l o oo o o 

is chosen to minimize |D(f)| or [j*?(D(f))] where 

M (f , ) = (1 - a ) M(f ) + a M(co ), 0 < a < 1 (3.41) 

1 o o o o o 

It can be shown that such an a exists . 

o 

(5) Repeat steps (2) - (4) until the desired accuracy is obtained. 


28 























Remarks 


1 . If k is the number of frequencies in the optimal design , then 


ml < k < -On + 1) 


(3.42) 


where 



m 


is the smallest integer higher than ^ 


2. The function \|/ has many local maxima. It is computationally time 
consuming to find C0 q where xj/ (co , f) is maximum. In the computer 
implementation of the algorithm, we consider finite but large N and 
search through all values of \J/(n) to find the maximum. Most stable 
systems of interest are low pass filters . Thus, in most cases, it is 
possible to find a subset of [-7C, 7i], where the search need be carried 
out . 


3. The termination criteria are the same as in Section 3.3.2 for the time 
domain case . 


Practical Considerations in the Computation of Optimal Input 


[27] 


The above algorithm will produce an optimal input design with a sufficient 
number of iterations. However, at each iteration, the procedure adds one point 
to the spectrum of the input. There are many inputs with unit power leading to 
the same information matrix. From a practical point of view, it is desirable to 
have as few frequencies in the optimal input as possible. The fewer the frequen- 
cies , the easier it is for the aircraft pilot or input generating system to reproduce 
the desired input. During the computation, a few steps can be taken to reduce the 
number of points in the spectrum . Suppose the normalized input at any stage has 
k frequencies ux with power (i = 1,2, . . . ,k) . 

(1) All frequencies which are "close" to each other can be lumped into 
one frequency . We consider two frequencies close to each other if 
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they are less than A 00 apart. Suppose q frequencies 00. are within 

1 * 

a band Aco wide. Then they can be replaced by one frequency CD 
* 

with power a where 
* *1 * 

a = S a. (3.43) 

i=i 1 

and 

* 1 ^ * 

co = — I a. oo. (3.44) 

a i=l 1 1 

(2) From this new normalized input, all frequencies ooj with power less 
than a threshold a' are chosen. These frequencies are not within 
Aco of any other frequency in the normalized input. They are 
dropped. The remaining frequencies do not form a normalized de- 
sign, so the design is renormalized. 

Steps (1) and (2) are carried out ensuring that the design does not become 
degenerate. This "practicalization" requires judgment of Aoo and a. 

3.6 OTHER CONSIDERATIONS IN THE DESIGN OF OPTIMAL INPUTS 

The stability and control derivatives determined from aircraft parameter 
identification are used for several purposes, for example, handling qualities 
specification, control system design, and simulator aerodynamic coefficient values . 
The ultimate objective of parameter identification enters into the considerations 
for the choice of optimal inputs to identify unknown system parameters . 

3.6.1 Primary and Secondary Parameters 

In many situations of practical interest, it is desirable to obtain accurate 
estimates of only a subset of all unknown parameters . Let the first k parameters 
be of primary interest and the remaining m-k parameters be of secondary 
interest. There is no direct incentive to obtain good estimates of the last m-k 
parameters . Inputs which give outputs sensitive to primary parameters should 
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be used. The secondary parameters are estimated only to the extent that they help 
reduce uncertainty in primary parameters resulting from errors in secondary 
parameters . 

This consideration is simple to handle when the optimality criterion is a 
weighted trace of the dispersion matrix. The weighting matrix should be chosen 
as 


m-k 


In the determinant criterion by selecting the function <p, 


cp = Tr(W M) 
T ps 


W 


D kk 

D kfi 

1 

o 



, k + fi = m 


-1 


(3.45) 


(3.46) 


(3.47) 


the is minimized. In other words, the input is such that it minimizes the 

area of the cross section of uncertainty ellipsoid mode by a k dimensional hyper- 
plane in the space of primary parameters . 


3.6.2 Technique for Evaluating Parameter Identifying Inputs 


Various techniques can be used to determine how effective an input is in 
identifying system parameters . The accuracy with which the parameters can be 
estimated depends not only on the inputs but also on the data reduction technique 
used . 


The most definitive way of determining the usefulness of an input is to use 
it with the system for which it is designed, and then to compare the resulting 
parameter estimates with the ones obtained using other inputs . If this is not 
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possible, simulation data must be used. There are two more basic methods to 
evaluate the simulation data: 

(1) Monte Carlo Simulation . This method consists of generating typical 
measurement time histories using simulations of system equations . 
The parameters are determined from each of these time histories by 
a minimum variance estimator. The mean and variance of the esti- 
mated parameters can be determined. 

/ 

(2) Information Matrix Method . It is known that the diagonal terms on 
the inverse of the information matrix are the covariance of the para- 
meter estimates if an efficient estimator is used. Information matrix 
can predict the parameter estimation errors without actually generat- 
ing a number of time histories . 

Monte Carlo simulation is usually expensive because the parameters must 
be estimated for each simulation using a minimum variance estimator (this is 
necessary so that the effects of the data reduction method do not mask the identi- 
liability of a certain input) . The information matrix method can accomplish the 
same task much more quickly. Any difference in the covariance of parameter 
estimates predicted by the two methods is only due to the inefficiency of the data 
reduction method used in the Monte Carlo simulation. 
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IV . TIME DOMAIN SYNTHESIS OF OPTIMAL INPUTS 


4.1 INTRODUCTION 

The time domain computer program determines the optimal inputs' time 
history for a specified length of flight test. In general, the optimal input depends 
on planned flight test duration, available instruments and their accuracies, the 
best a priori estimates of unknown parameter values, and many other factors. 

The techniques, developed in Chapter III, can be used to design the optimal input 
based on a variety of criteria. 

The elevator deflection sequence, which gives good estimates of parameters 

[3] 

in the short period mode of a C-8 aircraft, is computed under a variety of cir- 
cumstances . These inputs are evaluated against conventional inputs and against 
each other, with simulated data. The accelerations and velocities are determined 
over a simulated flight test to ensure that the inputs are safe and implementable . 

Rudder inputs are designed to identify five parameters in the lateral dynam- 

[3] 

ics model of a Jet Star aircraft. Combined rudder and aileron inputs are also 

determined and compared with rudder inputs alone. 

The weighted trace of the dispersion matrix is considered as one of the best 
criteria and is used throughout this chapter. Its comparison with other criteria 
for practical examples is given. 

4 . 2 INPUTS FOR LONGITUDINAL SYSTEM 

The longitudinal perturbation motions of a Buffalo C-8 aircraft about a 
steady trim condition obey the following differential equations (units: meters, 
deg, sec) . 
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where 



(4.1) 


u is forward speed, 

0 is pitch angle, 
q is pitch rate , 
a is angle-of-attack , and 
6 is stabilator deflection. 


The poles of this system are at 

CO = -1.16 + . 62j short period mode 
(Dpk = - .0153 + .088j phugoid mode 

The first pair of complex roots corresponds to the fast, highly damped, short 
period mode and the second to the slow, weakly damped, phugoid mode. 

It is assumed that there are five unknown parameters , all in the short period 
mode (underlined) . It is well known that the two state (pitch rate and angle-of- 
attack) model of an aircraft is a good representation of the short period motion. 
This approximation is used to find the elevator deflection time history. The 
equations of motion become 

' -1.588 -.562 1 fql -1.66 ‘ 

+ 6 (4.2) 

1 -.737 j [a] .005 s 
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with poles at 


co = -1 . 16 + . 62i 
sp — J 


These are noisy measurements of q and a 



(4.3) 


The measurement noise is assumed white . The root mean square errors in the 
measurements of q and a are .70 deg sec ^ and 1.0 deg, respectively, and the 
sampling rate is 25 per second. 


A doublet input, shown in Figure 4.1, is used conventionally to identify 

these five parameters . Starting from this doublet, the input design program is 

2 

used to determine the optimal input for a 6 sec long experiment with 100 deg sec 
total input energy. The performance index is the trace of the dispersion matrix. 
The input at the end of each iteration step is shown in Figure 4.2. Fairly good con- 
vergence is obtained in three steps. Table 4.1 shows the standard deviation in 

* 

parameter estimates for each of these inputs . Also shown is the value of P which 
determines the component of new input added to the old input (see Chapter III) . 

P was taken to be one if the decrease in Tr(D) was more than 50%. It is clear from 
Table 4. 1 that the optimal input should give much better parameter estimates 
(e.g. , smaller standard deviations) than the conventional doublet input. 


* 

The standard deviation format of Table 4. 1 is used throughout this report. These 
quantities are the square roots of the diagonal elements of the Cramer-Rao lower 
bound (Equation (2.1) , p. 8) , and represent the lowest possible value of the 
parameter estimate standard deviations which can be attained using an unbiased 
and efficient parameter identification procedure. As explained in Section 3.7, 
these lowest bounds rather than parameter estimates based on individual runs are 
a meaningful comparison of different inputs , a better input giving a smaller lower 
bound. The nominal parameter values are given for reference, in terms of 
dimensional coefficients , C. (j = m , mg , etc.), to simplify the notation. 

i q ° e 


37 





Table 4.1. Standard Deviations of Parameter Estimates for Inputs at 

the End of Each Iteration 

Length of Simulated Flight Test = 6 sec 

2 

Total Input Energy = 100 deg sec 



Standard Deviations 

Tr(D) 

P 

c 

m 

q 

■ 


C 

m g 

s 

C Z 

o 

Iteration 0 
(Doublet) 

.219 

.362 

.326 

.0978 

.0957 

.304 


Iteration 1 

.294 

.0703 

.0529 

.193 

.0452 

.134 

1.0 

Iteration 2 

.113 

.0729 

.0595 

.0620 

.0421 

.0272 

.61 

Iteration 3 

.113 

.0676 

.0561 

.0672 

.0400 

.0264 

.43 

Parameter Value 

-1.588 

-.562 

-.737 

-1.66 

0.0 



[All in units of deg, sec] 


Figure 4.3 shows pitch rate, angle-of- attack, pitch acceleration, and verti- 
cal acceleration time histories for the optimal input and compares them to the 
corresponding time histories resulting from a doublet input. Peak values of accel- 
erations are higher for the doublet input. However, high accelerations last for a 
longer time when the optimal input is used. Also, the excursions in angle-of- 
attack are much higher for the optimal input. 

4.2.1 Prim ary/Secondary Cons ide r ati. ons 

An input is designed considering C as the primary parameter of interest. 

mq 

The variance on the estimate of C is weighted 100 times more heavily than the 

mq 

variance on and . The starting input is the optimal input 

a a 8 o 

s s 

when all parameters are equally important. The input obtained after one iteration 
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Change in a (deg) Pitch Rate (deg sec 



Figure 4.3a 


Comparison of Pitch Rates for Doublet and Optimal 
Control (Simulation) 



Figure 4.3b Comparison of Angle-of-Attack Variation for Doublet 
and Optimal Input (Simulation) 



Pitch Acceleration (deg sec ) Vertical Acceleration (meters sec 








is shown in Figure 4.4. It lobks similar to the input when all parameters are 
weighted equally . The standard deviations of parameter estimates for these two 

inputs are compared in Table 4.2. The standard deviations of C and C m de- 

q 8 

M s 

crease by about 3% while those on other parameters increase . This shows that to 
get a good estimate of C m it is necessary to have a good estimate of C m also. 


Also, C is an important parameter in the system. Even when all the parameters 


m 


are to be identified, almost the best possible estimate of is obtained. The 


considerations of primary and secondary parameters may be more useful when the 
two sets of parameters affect different modes . 


Using the new time domain method, the entire idea of primary/secondary 
derivatives could be recast in a more general framework where parameters are 
assigned different levels of importance through weights on their covariance. 



Figure 4.4 Optimal Elevator Deflection (C Only Primary Derivative) 

m q 
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Table 4.2 Comparison of Standard Deviations on Parameter Estimates With 
All Parameters Equally Important vs . C Primary 

m q 

Length of Simulated Flight Test = 6 sec 

2 

Total Input Energy = 100 deg sec 



Standard Deviations 

Tr (D) 

c 

m 

q 

C 

m a 

C z 0 

C 

m § 

s 

C Z 

o 

s 

Parameters 

Equally 

Important 

.113 

.0676 

.0561 

.0672 

.0400 

.0264 

C Primary 

m q 

. 110 

.0747 

.0615 

.0651 

.0415 

.0275 

Parameter 

Value 

-1.588 

- . 562 

-.737 

-1.66 

0.0 



[All in units of deg , sec ] 


4.2.2 Comparison to Inputs Based on Maximizing the Trace of the Information 
Matrix 

Before the computer program was written to minimize the trace of the dis- 
persion matrix, many inputs were designed using the trace of the information 

matrix as the optimality criterion. Three such inputs 4, 6, and 12 sec long, 

2 

each with total input energy of 100 deg sec, are shown in Figure 4.5. Table 4.3 
compares these three inputs with each other and with the 6 sec. long time domain 
input obtained based on the Tr(D) criterion. It is apparent that the trace of the 
dispersion matrix is a much superior criterion. 

4.3 INPUTS FOR THE LATERAL SYSTEM 

The equations of motion for lateral motions of one version of a Jet Star flying 

[251 

at 573.7 meters/sec. at 6,096 meters are 1 (all in units of deg, sec ) 
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Figure 4.5 Longitudinal Inputs for Different Durations of Flight 

Test to Maximize Tr(M) 



where p is sideslip angle, cp is roll angle, and p and r are roll and yaw rates, 
respectively. Aileron and rudder are two control inputs. These are noisy 
measurements of the four states 
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(4.5) 


The noise in measurements is white and Gaussian with root-mean-square values 


rip 1 deg 

Tltp .5 deg 

_! (4.6) 

r|p .71 deg sec 

r)^ . 71 deg sec ^ 

and the sampling rate is 25 per second. The poles of this system are at 
-.0511 + 1.78 j Dutch -roll 

-1.12 Roll (4.7) 

-.00667 Spiral 

The inputs are designed to identify the parameters which predominantly affect 
the Dutch-roll mode (underlined) . 


4,3.1 Rudder and Aileron Inputs When All Parameters Are Equally Important 

The optimal rudder input to minimize the sum of dispersions of these five 

parameters is determined and is shown in Figure 4.6. The duration of the simu- 

2 

lated test is 8 sec. and the input energy is 100 deg sec. Figure 4.7 shows a 
conventional doublet input with equal energy . The comparison of standard devia- 
tion on parameters resulting from these inputs is given in Table 4.4. The optimal 
input results in better parameter estimates than the doublet, based on comparing 
the standard deviations . 
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Rudder Deflection (deg) , Rudder Deflection (deg) 
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[All in units of deg, sec. 



The state and lateral acceleration simulated time histories are compared in 
Figure 4.8. The optimal input results in large excursions in states and accelera- 
tions . This is because the Dutch-roll mode has low damping and optimal input 
continuously excites the responses. 

New simultaneous rudder and aileron inputs are designed to identify these 

2 

parameters . The inputs with combined energy of 100 deg sec are shown in 
Figure 4.9. The aileron input amplitude is very small. The estimates resulting 
from this simultaneous input are presented in Table 4.4. There is a very small 
improvement over the single rudder input case, as would be expected since the 
rudder input is much more effective in estimating these five parameters than the 
aileron input. Larger aileron deflections can be obtained by placing separate 
energy constraints on the aileron and the rudder . 

There are eight linearly independent columns in the controllability matrix 
of (Fq, Gq). Therefore, eight linear differential equations must be solved to 
obtain all sensitivity functions (see Section 3.3.3 and Appendix D). However, the 
last column is almost linearly dependent on the remaining seven. The optimal input 
is designed through propagating just seven equations . The standard deviations of 
parameters resulting from this input are given in Table 4.4. There is an insig- 
nificant difference in this input and the input based on using all eight sensitivity 
equations. This approximation brings about considerable saving in computation 
time. It could prove extremely useful in computation of optimal inputs for practical 
high order systems . 


4.3.2 Primary/Secondary Derivatives 


A rudder deflection sequence is determined to primarily identify C 


C and C 
n n c 

r 8 


and consider C 


n p 

and Cg as unknown secondary coefficients. As 


y p ~p 

X 

in the previous case, this is accomplished by weighting the error covariance on 
primary parameters 100 times as on secondary parameters . The input is shown in 
Figure 4. 10 (one iteration starting from the optimal input of Figure 4.7) . The inputs 

and C 


considering all parameters primary and considering only C , C 

n r n p 


n c 
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Figure 4.8b Lateral Acceleration for Doublet and Optimal 

Input (Simulation) 








Yaw Rate (deg sec ) era Roll Rate (deg sec 



Figure 4.8d Yaw Rate for the Optimal Rudder Input (Simulation) 

















primary are quite similar. The parameter error covariances are compared in 
Table 4.5. There is only a moderate improvement in the accuracy of the primary 
parameters . 


4.4 CONCLUSIONS 

Inputs are designed for both longitudinal and lateral systems based on a 
variety of criterion functions. For these inputs, the following conclusions are 
reached: 

1. Inputs based on the trace of the dispersion matrix criterion give 
lower derivative estimate error standard deviations than those based 
on the trace of the information matrix criterion. 

2 . Low order approximations to system representations can usually be 
made in the computation of the optimal input (two state representation 
of short period mode gives almost as low estimate error variances as 
the more complete, four state longitudinal model). 

3. For weakly damped systems, it may be necessary to place direct state 
constraints to avoid excessive state excursions. It is noted, however, 
that weakly damped modes causing large aircraft motions are accept- 
able as long as the pilot can tolerate and control them. 

4. The inputs designed to identify primary parameters alone did not show 
much improvement in the estimation accuracy of these parameters 
ove;r the inputs designed to identify all parameters . A more careful 
separation of primary and secondary parameters may be required in 
the design of optimal input to gain maximum benefits from these 
considerations. In the procedure used here, it is assumed that there 
is no knowledge about the secondary parameters before the flight 
test. Further work is required to include cases in which partial in- 
formation about the secondary parameters can be obtained from an 
independent source . 
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Table 4.5 Comparison of Standard Deviations on Parameter Estimates for Inputs 
(All Parameters Equally Important vs . C , C , and C Primary) 
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[All in units of deg, sec 
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V . FREQUENCY DOMAIN SYNTHESIS OF OPTIMAL INPUTS 


5 . 1 INTRODUCTION 

The frequency domain inputs are designed on the assumption of steady state 
and are "optimal 11 only when the flight test duration is long as compared to the 
time constants of modes of interest. In aircraft applications, this may not be true. 
Nevertheless , frequency domain methods play an important role in the design of 
optimal inputs even for "short" experiments. The input spectrum, which specifies 
steady state optimal input, can be computed quickly. Subsequently, the time trace 
obtained by adding different frequency components with a random phase relation- 
ship is usually a good first pass at the optimal input time history for short experi- 
ments. It can be used as such if the flight test is not too short and is an excellent 
starting input for the time domain program. 

A detailed study is made of the effect of short data length on the optimality 
of frequency domain inputs . Input sequences based on different experiment dura- 
tions are generated from the spectrum and evaluated on the system, which starts 
from zero initial condition and does not reach steady state. 

r?7l 

5 . 2 FREQUENCY DOMAIN OPTIMAL INPUTS L J 

Under the present effort, a computer program is written to design optimal 
inputs in frequency domain for single input systems . It is used to find optimal 
elevator and rudder input spectra to estimate parameters in longitudinal and 
lateral modes, respectively. 

5.2.1 Longitudinal System 

As a first example, an optimum elevator deflection spectrum is derived to 
identify parameters in the longitudinal short period mode of a C-8 aircraft. The 
state and measurement equations are given in Section 4.2. Two criteria of opti- 
mality are tried. They are: (a) min |D| , and (b) min Tr(D). To use the computer 
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program, the discretization step is chosen as .02 sec and the number of points is 


4000. Since there are five parameters and two outputs, is two. Thus, the 

minimum number of frequencies for a nondegenerate design is two and the maximum 

5x6 


number of frequencies required in an optimal design is 
(Equation (3.42)). 


, i . e . , fifteen 


To minimize the determinant of the dispersion matrix, the initial input is 
selected to have two frequencies at 0 cps and at .125 cps with equal power. Fig- 
ure 5.1 shows the spectrum of the elevator deflection input after each iteration. 
Notice that during some iteration steps, the program puts more power at already 
chosen frequencies . After eight iterations, the change in the determinant of the 
dispersion matrix is less than . 1% from the previous step. There is a total of 
seven frequencies in the final input spectrum. 

This input has interesting characteristics . The spectrum is divided into 

two parts: A low frequency input to identify gains and a high frequency input 

to identify natural frequency and damping, etc. The higher frequency input 

occurs around the natural frequency, which is reasonable. Table 5.1 shows the 

standard deviations (lower bound) on parameter estimates for an average 12 sec 

long experiment when the system is in steady state and the total energy in u dur- 

2 

ing this time is 200 deg sec. Also shown are the trace and determinant of the 
information matrix and the trace of the dispersion matrix. Next, the frequencies 
close to each other are lumped and ones with too little power are dropped to get 
the spectrum of Figure 5.2. The standard deviations on parameter estimates for 
this simplified design in steady state for the same experiment duration and input 
energy are also shown in Table 5.1. There has been an improvement in the deter- 
minant of M and trace of D, showing the value of the simplification. 

The frequency domain inputs are designed with the assumption of steady 
state. If the flight testing time is short, the aircraft does not reach a steady 
state. To find the true information matrix for a 12 sec long test starting from 
zero initial conditions, a time domain input based on the frequency spectrum and 
same average power is generated and is shown in Figure 5.3. This time domain 
input is not unique since the initial phase relationship between the sinusoidal 
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Steps in the Computation of Optimal Input (min |D|) to Identify 
Parameters in the Short Period Mode of a C-8 Aircraft 
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[All in units of deg, sec 
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Figure 5 .3 Elevator Deflection Time History Based on Spectrum of 

Figure 5.2 (min |D|) 








waves is arbitrary. Table 5.1 gives the parameter standard deviations and trace 
and determinant of information and dispersion matrices when initial phases are 
chosen at random. The parameter standard deviations deteriorate by 5% to 15%. 

A better result could be obtained by optimizing the initial phases . 

The basic and simplified elevator deflection spectra for the Tr(D) criterion 
are shown in Figure 5.4. The frequencies in this input are in the same range as 
in the input for the |D| criterion. Table 5.2 shows standard deviations in para- 
meter estimates for an average 12 sec period in the steady state when the total 

2 

input energy is 200 deg sec. The simplified design is a little poorer than the 
basic design. Some parameters have a higher standard deviation than in the |D| 
case, while others have lower standard deviations. Again, a 12 sec long time 
domain input (shown in Figure 5.5) is generated and the standard deviations for 
this input with the system starting from zero initial condition are given in Table 
5.2. There is only a moderate deterioration from steady state value. 
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Figure 5.5 Elevator Deflection Sequence to Minimize Tr(D) 
«' 1 


The root-mean-square (RMS) state deviations are computed for the frequency 

2 

domain inputs . For an average input power of 16.67 deg , the RMS states for the 
inputs of Figures 5.2 and 5.4 are shown in Table 5.3. The values look reason- 
able . 


Table 5.3 RMS State Deviations for Frequency Domain Inputs 

(Simulation) 


Criterion 

Pitch Rate 

Angle-of- Attack 

|D| 

3 . 19 deg sec ^ 

2.95 deg 

Tr(D) 

3.05 deg sec ^ 

2.52 deg 


64 





5.2.2 Effect of Short Data Length on Performance of Frequency Domain Inputs 


If the system starts from zero initial conditions, the performance of an input 
of finite length is poorer than predicted on the assumption of steady state. To 
determine the duration of the experiment when this approximation becomes serious, 
time traces of elevator deflection 4 to 12 sec long are obtained based on the simpli- 
fied spectrum of the Tr(D) criterion (Figure 5.4) . Each of these inputs is used 
with the state and measurement equations and the resulting information and disper- 
sion matrices are determined. Table 5.4 gives the ratio of standard deviations on 
parameter estimates for these finite inputs (with the system starting from zero 
initial conditions) to the standard deviations in steady state (for the same average 
input power and experiment duration) . Trace and determinant of the information 
and dispersion matrices are also compared using simulated data. The asymptotic 
value of these ratios for long experiments is one. For experiments shorter than 
8 sec, the deterioration is serious. The inputs are good for experiments longer 
than 10 sec. This corresponds to about two cycles of the natural short period 
mode . 

5.2.3 Comparison With Conventional and Optimal Time Domain Inputs 

The inputs, which minimize the sum of variances of five parameters in the 
short period mode of a C-8 aircraft for a 6 sec long experiment are given in 
Chapter IV. There, the conventional doublet input is also given. 

Table 5.5 shows a comparison of standard deviations on parameter esti- 
mates for the frequency domain input, the optimal time domain input, and the 

doublet. The steady state frequency domain value is a lower limit on Tr(D) for 

2 

an input with 100 deg sec input energy in a 6 sec long flight test. The time 
domain input is optimized for a 6 sec long experiment and gives a much better 
result than the time trace from the frequency domain input. Nevertheless, the 
input resulting from the frequency domain approach is superior to a doublet. As 
mentioned before, this would be an excellent first pass at the optimal input and 
is useful for starting the time domain program. 
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Table 5.5 Comparison of Time Domain and Frequency Domain Approach 

(Tr(D) Criterion) 
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Figures 5.6 and 5.7 show the simulated pitch rate and angle-of-attack 
variation during the 6 sec long experiment for both inputs . The peak values are 
comparable . 

5.2.4 Frequency Domain Lateral Inputs 

The lateral system of Section 4.3 with unknown C , C n , C , C and 

yp g p n p n r 

C is used as the example to determine the optimal rudder input spectrum to mini 

n 8 

r 

mize the trace of the dispersion matrix. The optimum input spectrum has two 
frequencies: at 0 cps with 12% of total input power and at .285 cps with 88% of 
input power. The second frequency is very close to the natural frequency of the 
Dutch-roll mode. Since the Dutch-roll mode for this aircraft has low damping, this 
input would produce large state excursions in steady state. This occurs because 
there are no constraints on the state variables . 

Table 5.6 shows a comparison of standard deviations for this frequency 
domain input and the time domain input designed in Chapter IV for an 8 sec long 
simulated flight test. Because of low damping, the system is far from steady state 
for the duration of the experiment. Standard deviations on parameter estimates 
predicted on the assumption of steady state are too optimistic. 

5,3 CONCLUSIONS 

Though the frequency domain inputs are optimal only if the flight test is 
long, they are useful even for short experiments where steady state conditions 
cannot be reached. These inputs could be used as they are for flight tests which 
are longer than two cycle times of the mode of interest. Since frequency domain 
inputs can be obtained in a much shorter computation time , this technique may 
have advantages in real-time and on-line applications. These inputs are excellent 
for starting other more complicated algorithms, for example, the time domain 
input design. 
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Figure 5.6 Pitch Rate Time Histories for Time Domain and Frequency 
Domain Optimal Inputs (Simulation) 
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Figure 5.7 Angle- of- Attack Time Histories for Time Domain and 
Frequency Domain Optimal Inputs (Simulation) 
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limit on Tr(D); cannot be achieved starting from zero initial conditions. 



































A small input energy near the natural frequency of a lightly damped mode 
can cause large state excursions . Since the input design places a constraint on 
the input energy, there would be a concentration of power near the natural fre- 
quency, because this would increase the signal-to-noise ratio. If this is undesir- 
able, it can be avoided by putting state constraints or limiting the allowable fre- 
quencies in the input spectrum. 
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VI . EVALUATION OF TIME AND FREQUENCY DOMAIN OPTIMAL INPUTS 


6 . 1 BASIS OF EVALUATION 

The optimal inputs of Chapter III are evaluated for performance sensitivity 
when design conditions are not satisfied. There are two possible sources of error. 
The control surface deflection may be different from the optimum because of the 
pilot's inability to follow the input exactly. Secondly, the input design may be 
based on incorrect parameter values (based on incorrect a priori estimates) or 
overly simplified models . 

The optimal inputs are approximated by a series of four steps . A simulation 
of these step inputs is used on the system to evaluate parameter error covariances 
resulting from the measurements. This evaluates the degradation in performance 
of optimal inputs from errors in implementation. 

Next, it is assumed that the parameters of the system are 50% off from the 
design values. The optimal inputs based on design conditions are used under 
off-design conditions . A comparison of these inputs with conventional inputs 
determines the loss in input efficiency resulting from inaccurate knowledge of 
parameters . 

In Chapter IV, the optimal elevator inputs were computed based on the two 
state approximation of the short period mode in the longitudinal motions . Those 
inputs are used in this chapter on the full four state model to study the effects of 
model approximation on performance of the input resulting from it . 

6 . 2 LONGITUDINAL SYSTEM 

6.2.1 Approximation to Optimal Input 

Figure 6. 1 shows the optimal input of Chapter IV to identify five short 
period parameters of a C-8 Buffalo aircraft. This input is approximated by a 
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sequence of four steps with the same total energy. The approximated input is shown 
by the broken line in Figure 6.1. This input is used on the two state, short period 
approximation of the longitudinal equations of motion (Equation (4.2)) . Table 6.1 
compares parameter dispersions resulting from the optimal and the approximated 
(suboptimal) inputs . Some parameters have better estimates while others have 
poorer estimates . 


Table 6. 1 Comparison of Optimal and Approximated Inputs 

Duration of Simulated Flight Test = 6 sec. 

2 

Total Input Energy = 100 deg sec 



Standard Deviations 

Tr(D) 

c 

m 

q 

C 

m a 

% 

C 

m 6 

s 

c z 

s 

Optimal Input 
(Solid line , 
Figure 6 . 1) 

.113 

.0676 

.0561 

.0672 

.0400 

.0264 

Approximated Input 
(Broken line, 
Figure 6 . 1) 

.126 

.0590 

.0497 

.0807 

.0381 

.0299 

Parameter Value 

-1.588 

-.562 

-.737 

-1.66 

0.0 



[All in units of deg, sec.] 


The pitch rate and angle-of- attack for the optimal and the approximated 
input are compared in Figure 6.2. The maximum values are about the same in the 
two cases . 

6.2.2 Off-Design Parameter Values 

All five parameters in the short period model are increased by 50% of their 
initial values . This results in a system with a natural frequency of 1.86 rad/sec and 
a damping ratio of .94. It is more difficult to identify parameters of this system with 
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Figure 6.2b Angle-of- Attack for Optimal and Approximated 

Inputs (Simulation) 







the same input energy because of increased damping. Table 6.2 shows standard 
deviations of parameter estimates when the approximated input is used on this 
system with off-design parameter values . There is a sharp increase in estimation 
errors from design conditions, partly because of the higher damping ratio. Table 6.2 
also shows parameter error covariances with a doublet input. Next, the para- 
meters of the system are halved, reducing the natural frequency to . 76 rad/sec 
and the damping ratio to .77. The parameters in this system are easier to identify. 
The performances of the approximated input (broken line of Figure 6.1) and the 
doublet are given in Table 6.2. For the approximated input, the parameter error 
covariances are smaller than under design conditions. In both cases, the approxi- 
mated input compares favor ably with the conventional doublet input. Figure 6.3 
compares the simulated pitch rate and angle-of-attack for the approximated input 
under design and off-design conditions (all parameters halved) . 

6.2.3 Fourth Order Model 

The approximation to optimal inputs, obtained using the two state approxi- 
mation, is simulated on the four state longitudinal equations of motion (Equation 
(4.1)) . Again, there are measurements of only q and a. The measurement 
error and the sampling rate are the same as before. Table 6 . 3 is a comparison 

of the standard deviations on estimates of C , C , C , C- , and C n on the 

m m m e L 
q a 8 a 6 

^ s s 

assumption that the remaining parameters in the system are known. The estimates 
predicted by the four state model are better than the estimates predicted by the 
two state model. This is because there is an additional excitation of the short 
period mode from variations in forward speed. 

Simulated time histories of the deviations in forward speed, pitch rate, and 
angle-of-attack are shown in Figure 6.4. They are compared to the output of the 
two state model. There is an insignificant difference, except in forward speed 
which is assumed to remain constant in the two state short period approximation 
of the longitudinal equations of motion. 
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Buffalo 08 Aircraft 
Altitude = Sea Level 
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Figure 6.3a Pitch Rate Variation for Systems Under Design and Off- 
Design Conditions (Approximated Inputs) Based on Simulation 
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Figure 6.3b Angle-of- Attack Variation for Systems Under Design and 
Off-Design Conditions (Approximated Input) Based on Simulation 
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4 State 



2 State 



Buffalo C-8 Aircraft 
Altitude = Sea Level 
Velocity = 41.2 meters/ sec 
Input Energy = 100 (deg)^(sec) 


Time (sec) 


Buffalo C-8 Aircraft 

Altitude = Sea Level 

Velocity = 41.2 meters/sec 

2 

Input Energy = 100 (deg) (sec) 


2 State 


Time (sec) 


2 State 
4 

6 

Time (sec) 

\ 




h 4 State 

Buffalo C-8 Aircraft 


\ 

Altitude = Sea Level 
Velocity = 41.2 meters/ sec 
Input Energy = 100 (deg)^(sec) 


Figure 6.4 Comparison of State Excursions Predicted by Two State 
and Four State Models (Approximated Input) Based on Simulation 
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Table 6.3 Comparison of Standard Deviations on Parameter Estimates 
Predicted by Two State and Four State Models 

Length of Simulation = 6 sec 

2 

Total Input Energy = 100 deg sec 


Model 

Standard Deviations 

Tr(D) 

c 

m 

q 

C 

m a 

% 

C 

m 5 

s 

c z 

o 

s 

Four State 
[Equation (4.1)] 

.115 

.0560 

.0499 

.0752 

.0378 

.0260 

Two State 
[Equation (4.2)] 

.113 

.0676 

.0561 

.0672 

.0400 

.0264 

Parameter Value 

-1.588 

-.562 

-.737 

-1.66 

0.0 



[All in units of deg, sec ] 


6.3 LATERAL MODES 

The lateral system equations and the optimal input to obtain good estimates 
of C , C 0 , C , C , and C are presented in Section 4.2. As for the lateral 

y p 8 p y k \ 

system, this input is approximated by a series of four steps and is shown in thq, 
broken line in Figure 6.5. Table 6.4 compares the parameter error dispersions 
for the optimal and the approximated inputs. There is less than a 20% increase in 
the trace of the dispersion matrix (this corresponds to an average increase of 10% 
in standard deviation on parameter estimates) . 

The parameters are next increased by 50%. This increases the damping 
and natural frequency of the Dutch-roll mode significantly. Table 6.5 compares 
the standard deviations on estimates with optimal input and doublet under off- 
design conditions. Though the standard deviation on parameter estimates rises 
considerably, the optimal input is still much better than a doublet. It is important 
to note that a part of the degradation in estimation accuracy is due to increased 
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Table 6.4 Comparison of Optimal and Approximated Rudder Inputs 


Duration of Simulated Flight Test = 8 sec 

2 

Total Input Energy = 100 deg sec 


Input 

Standard Deviations 

Tr(D) 

C 

Y P 

C fi 

P 

C 

n P 

c 

n 

r 

C 

n 8 

r 

Optimal 

.00880 

.0204 

.00277 

.00880 

.00860 

.000648 

Approximated 

.00933 

.0221 

.00310 

.00977 

.00908 

.000764 

Par ameter 
Value 

-.119 

-4.43 

2.99 

-.178 

1.55 



[All in units of deg , sec ] 


damping ratio under off-design conditions which make parameters less identifiable 
from a given input energy . 

6.4 CONCLUSIONS 

Large errors in the implementation of optimal inputs result in small deter- 
ioration in parameter estimation accuracy. In particular, it is possible to approxi- 
mate the optimal inputs by a series of steps at a small cost in increased standard 
deviations on parameter estimates. This is true both of a system with low damping 
and of a system with high damping in the mode of interest . 

Changes in parameter values from design conditions investigated here 
changed the damping and natural frequency of the mode of interest considerably. 
For the same input energy, the parameter of the resulting system can be estimated 
with either higher or lower accuracy, depending mainly on the increase or de- 
crease in damping and to some extent on changes in natural frequency . The opti- 
mal input is sensitive to off-design parameters; however, with a 50% change in all 
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parameters from their design values, the optimal inputs are at least as good as 
conventional doublet inputs . The deterioration in performance of the optimal 
input is more severe when the damping is low with a priori estimates of parameter 
values but is higher for true parameter values . 
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VII . SUMMARY AND CONCLUSIONS 


A method for designing inputs for linear time-invariant systems is developed 
for application to the identification of aircraft stability and control coefficients. 

The method consists of two complementary techniques: 

1. The time domain technique can be used for any linear dynamic 
system under various conditions and all durations of flight test time . 

It is general, but computationally complicated. 

2. The frequency domain technique gives an input spectrum, the time 
traces of which are good inputs for parameter identification for flight 
tests which have long allowable data lengths compared to the time 
constants of the mode of interest. This method could also be used 

to determine a set of specific frequencies where the flight test may 
be carried out, obviating the need for lengthy in-flight frequency 
sweeps . 

Both of these techniques can be used for input designs when test constraints 
allow for relatively long time aircraft response to such inputs . The time domain 
method gives the best results and should be used whenever possible in spite of 
its computational complexity, because computer time is usually a less important 
constraint than flight test time. The frequency domain method has application in 
obtaining fairly good inputs quickly for long experiments, and in finding starting 
inputs for the time domain program . 

As the limitation on allowable aircraft response time becomes stricter, the 
time domain method is the desirable approach . For example, a test at an extreme 
flight condition (e.g., high angle- of- attack, sonic transition, etc.) may preclude 
extensive times for response data acquisition. At such points, the time domain 
capabilities are best used to generate shorter inputs which supply maximum 
parameter estimate accuracy. It is also noted that considerations of cost and/or 
fuel restrictions may be achieved by using efficient inputs requiring less data 
acquisition time . 


87 



The improved applicability of these optimal input design techniques to flight 
test is achieved because of the resolution of several computational problems 
which have previously hindered such application. The solution to these problems 
has been attained by the following accomplishments: 

(a) Design of inputs based on a wider class of optimizing criteria 
which address more directly the requirements for which the para- 
meters are identified. The criteria are based on the variances of 
parameter estimates and include determinant and trace of the disper- 
sion matrix. Other criteria may also be used (cf . , Appendix A 

and Section 3 . 4 . 2) . 

(b) Increase in the number of parameters to be identified by design 

of inputs for a larger number of parameters than previously computa- 
tionally feasible. This benefit is obtained by implementation of two 
new algorithms . These are: 

(i) Sensitivity Function Reduction : The technique of sensitivity 

function reduction lowers the order of the computational. prob- 
lem which must be solved in the time domain input design 
method. For example , in a system with 4 states, 2 inputs, 

4 measurements, and 15 unknown parameters, the number of 
differential equations is reduced from 64 to a majqmum of 12. 
Since the computation time varies as the cube of the order of 
the problem, this sensitivity reduction method would reduce 
computation time by a factor of more than 100 (Appendix D 
and Section 3.4.3) . By using sensitivity functions reduction, 
it has become numerically feasible to compute inputs for 
high order systems with many unknown parameters ■ 

(ii) Eigenvalue-Eigenvector Decomposition Method to Solve the 

Two Point Boundary Value Problem : Properties of the transi- 

tion matrix occurring in the time domain method are used to 
develop an eigenvalue-eigenvector decomposition solution 
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technique. This has been found to be quite efficient and gen- 
eral (Appendix C and Section 3 . 4 . 4) . 

(c) Frequency domain input design specification is desired specifically 
for tests where sinusoidal input generators are available. By making 
the steady state assumption, the input design procedure can be 
simplified considerably. In this method, the frequency domain repre- 
sentation of a system is used to determine the optimal input spectrum. 
This technique is useful for designing long inputs (Appendix B and 
Section 3.5) and obtaining first estimates of inputs for short experi- 
ments . 

(d) Specially designed inputs for primary /secondary derivative con- 
sideration allow flexibility in tailoring inputs to identify specific 
important parameters. Alternate formulations where a varying amount 
of confidence in "known" parameters exists can be easily incorpor- 
ated into the methods (cf.. Section 3.6,1). 

The evaluation and verification of these benefits have been performed by 
the design of optimal longitudinal and lateral inputs and subsequent testing on 
simulated data. These evaluations have led to the following conclusions: 

(a) Inputs based on the dispersion matrix criteria Tr(D) and |D| give 

lower parameter error covariances than inputs based on the informa- 
tion matrix criterion Tr(M) (Section 4.2.4). For example, a longi- 
tudinal input based on Tr(D) gave C and C estimate error 
^ m m 0 

q 8 

n s 

deviations which were approximately one-third of that for the input 
given by the Tr(M) criterion (cf., Table 4.3). This conclusion is 
significant because it confirms the advantages which were anticipated 
by reformulating the problem to be more applicable to flight test 
requirements . This reformulation involved using the dispersion 
matrix directly, a previously intractable computational procedure 
for a large number of parameters . 
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(b) 


Optimal inputs give lower parameter error covariances than con- 
ventional control doublets (for the same input energy) . The doublet 
was chosen as the representative input because of its wide use as a 
basic flight test maneuver. In all cases , the optimal input gave 
derivative estimate accuracies almost 50% or better than the doublet 
(c.f. , Table 4. 1 for longitudinal and Table 6.5 for lateral examples) . 
The greatest improvement in accuracy through optimal inputs is 
obtained in derivatives affecting highly damped modes , and this is 
one singular benefit of such optimal inputs (cf. , Section 4.2 and 
4.3) . Even under off-design conditions, where the assumed deriva- 
tives for calculating the optimal input were as much as 50% off from 
the "true" values, the doublet gave higher error covariances (cf., 
Sections 6 . 2 and 6.3). 

(c) Optimal inputs may often be approximated by a series of steps 

without significant increase in derivative estimation errors (Sections 
6.2 and 6.3). For example, the increase in the sum of parameter 
error covariances was less than 15% when the optimal input was 
approximated by a series of four steps for both longitudinal and 
lateral simulated responses . 

Several useful guidelines have been established for use of the time and 
frequency domain input design techniques . Such guidelines are useful for reduc- 
ing the computation time required for designing optimal inputs . These include: 

(a) Based on the computed results, it is suggested that frequency domain 
inputs be designed when the experiment duration is longer than 
two cycle times of the mode of interest (Section 5.2.2) . 

(b) There are several practical approximations possible in the computa- 
tion of optimal inputs . These simplifications reduce design time with 
only a slight loss in accuracy. 
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(i) 


Inputs can be designed based on simplified models . A two 
state representation of the short period mode is adequate in the 
design of inputs to identify the five parameters which affect 
this mode (Section 6.2.3). 

(ii) Sensitivity equations can, sometimes, be reduced below the 

minimum number required to obtain the exact value of all sensi- 
tivity functions . This approximation was made for the lateral 
system in Section 4.3 and compared to the exact case in 
Table 4.6. 

(c) Direct state constraints in the design may be desired , particularly 
in weakly damped modes to avoid excessive aircraft responses 
(Section 4.3) . As long as the pilot can control the aircraft, however, 
such motions will not limit the use of the inputs studied in this report. 

These results are preliminary, but do indicate that the approaches developed 
are suitable for flight test applications. 
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APPENDIX A 

INPUT DESIGN IN TIME DOMAIN 


A. 1 PROBLEM STATEMENT 

Consider a continuous time varying system 

x = F(t, 0)x + G(t, 0)u 
x(0) =0 0 < t < T 

where 


(A . 1) 


x is a n x 1 state vector 
u is a q x 1 control vector, and 

F(t, 0) and G(t, 0) are n x n and n x q matrices, which depend on m 
unknown parameter 0. 


There are noisy measurements of some linear combinations of state 

y = H(t, 0)x +v z + v (A. 2) 

y is a p x 1 output vector and v is p x 1 white noise vector with zero mean and 
known covariance R(t). 

The problem is to choose u from a class of second order processes inde- 
pendent of v to obtain good estimates of parameters from the measurements of 
u and y. The class of inputs is such that 


u (t) u ( t ) = 1 


(A. 3) 
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A . 2 INFORMATION MATRIX 

The information matrix for parameters 0 is 

M= / (lf) T If 0‘ 

*o 

(A. 1) can be solved for state vector x. 


x(t) = / <Kt, t) G(t, 0) u(T) dx 

o 


where <p(t, t) is the transition matrix and follows the differential equation 


= F(t ’ 0) c * )(t ’ T) 


<J>(x, t) = I 


Vector z can be written as 


z(t) = / H(t, 0) <j>(t, x) G(t, 0) u(x)dx 
o 

q r 

= 2 / H(t, 

i=l J_ 


Therefore , 


3z(t) 

30 


q 

E 

i=l 


0) <|>(t, x) G^x, 0) u ± (x)dx 


j Jq (H(t, 0) <Kt, x) G i (x, 0)} u i (x)dx 


(A. 4) 


(A. 5) 


(A. 6) 


(A. 7) 
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/ t 

A. (t , t) u.(x)dx 
1 


(A. 8) 


Using (A. 7) in (A. 4), 
.T 


f , f t 
u - Jo di K A * (t - T> 


(T)dT R 1 (t) 


q 

z 

j=l 


A (t, s) u. (s)ds dt 
J J 


n r q r T 

Z Z / |A T (t, t) R \t) A.(t, s dt/u.(x)u. 

i=l j=l J 11 J 1 J 


(s) dx ds 


o o 


sup (x , s) 


where sup(x, s) is the higher of X and s. This can be written as 


M = 



q 

Z M. . (x, s) u.(x)u.(s)dx ds 
i,j=l 1J . 1 3 


where 


s) = 1/2 J ^Ct, T) R 1 (t) A.(t, s) 


sup (x , s) 


(A. 9) 


+ aJ ( t, x) R~ 1 (t) A.(t, s)}dt 


It is clear that 

M ij (T » s) = M ji( x » s ) = m£.(s, x) 


(A. 10) 


(A. 11) 
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The next sections show how one can work with the three criteria of optimal- 
ity outlined in Chapter in . 

A . 3 MAXIMIZATION OF A LINEAR FUNCTION OF THE INFORMATION MATRIX 


From (A. 9), a linear functional SB (see Section 3.3 for definition of SB) of 
the information matrix is 


J 1 A SB(M) 



SB(M.. (t, s)) u.(t) u. (s) dx ds 
13 i J 


where 



(x) P(t, s) u (s) dx ds 


(A. 12) 


P ±j< T » s ) s)) 


(A. 13) 


P(t, s) is a symmetric positive definite matrix. To maximize a linear function of 
the information matrix for the given class of inputs, it is necessary to maximize 
(A. 12) under the constraint (A. 3) . It is straightforward to show that this is 
achieved by solving the following eigenvalue problem: 



o 


P(t, s) u(s) ds = Au(t) 
0 < t < T 


(A. 14) 


The maximum eigenvalue X is the maximum value of the trace of the information 
matrix and the corresponding eigenvector is the desired optimal input. Mehra 
shows this eigenvalue problem can be recast into a linear two-point boundary 
value problem. The solution technique to this problem is given in Chapter III. 


[ 22 ] 
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A . 4 DETERMINANT OF THE DISPERSION MATRIX 


In this and the next sections, we prove some important theorems. These 
lead to a computation algorithm presented in Section 3.4.3. 

Theorem A . 1 : A necessary condition that an input u (t) minimizes |D*| is that 

* * 

(i) u is an eigenvector of P with eigenvalue m. 

k ^ 

(ii) Any other eigenvalue, with eigenvector u , of P follows the 

inequality 


A k < m + 2 Tr{(M X M ) } 


(A. 15) 


where 


p ij (T> s) = Tr(M * 1 M ij (T ’ s)) 


(A. 16) 


M 


*k 


r T rT 

= I 2 M (x, s) u.(t) u^(s) dx ds ( A - 17 ) 

o •'o i ’ j=1 J 1 3 


and 


P (x, s) u^s) ds = A. it<T) 


(A. 18) 


Proof: (a) From (A. 11) 


* 

M 



q 

£ M..(x, s) u*(x) u*(s) dx ds 

i,j=l 1 J 
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Tr (I) 



Q * — 1 

Z Tr(M s)) u*(x) u*(s) dx ds 

i,j=l 1 J 


m = 



u* T (x) Ptr, 


s) u\s) dx ds 


(A. 19) 


Therefore, X >m (A. 20) 

max — 

* * 

Equality holds if u (t) is an eigenvector of P corresponding to 
its maximum eigenvalue m. 

5}C 

(b) Consider an input u (t) which minimizes |D| . Since logarithm is a 

* 

monotonic function, u (t) must also minimize log |D| . Let V be 

* k K 

an eigenvalue of P and u (t) the corresponding eigenvector . 

Consider an input u, 

u(t) = a u*(t) + 8 u k (t) (A. 21) 

The energy constraint (A. 3) requires 


a 2 + 


8 + 2a8y = l 


where 


Y 


-f 


u (t)u (t)dt 


giving 


9a _ (8 + ay) 

98 (a + By) 


(A. 22) 


(A. 23) 


(A. 24) 


The information matrix for input u(t) is 



(A. 25) 


2_ * 2 k *k 

M = a m + 3 M + 2a BM 


where 


■/’/' 
o o 


M k 


A 

. ? M ii^ T ’ U ^CT) u k (s) dx ds (A. 26) 
i»J=l J 1 J 


9 log |d 


&3 


3=0 


^ 1 ^ ^ 

= -Tr{M (-2yM + 2M k )} 


2y(m - X k ) 


(A. 27) 


Since u minimizes |D| and y can be positive or negative, for all 
* 

eigenvalues of P 


A, = m or y = 0 


(A. 28) 


5§C * 

If Y ^ 0 for one particular eigenvector of P (except possibly u ) , 

* 

it cannot be zero for any other eigenvector of P , since the eigen- 
vectors of a positive self adjoint kernel are orthogonal. Then all 
* 

eigenvalues of P must equal m. In the light of Equation (A. 19) , 
* * 

this implies that u is also an eigenvector of P with eigenvalue 
m . Then y = 0 for all other eigenvectors . 


If y = 0, then u is an eigenvector and the corresponding eigen- 

* * 

value is m. Thus, we show that u is an eigenvector of P with 
eigenvalue m. 


. . 

Since u minimizes |D| 
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= -Tr{-M* 1 (M* k M*~ 1 M* k + 2M* - 2M k ) } 


=4Tr{(M* V*) 2 } + 2m - 2 A, 


A 1 m + 2 Tr{ (M* V*) 2 } 


Theorem A, 2: A sufficient condition that an input u (t) is a stationary point for 

ii»* * 

|D| is that u be an eigenvector of P with eigenvalue m. 

Proof : Assume u is an eigenvector of P . Consider any input u (t) satisfying 

the energy constraint and define 


u(t) = au (t) + 3u m (t) 


Since the eigenvectors of P form a complete set and are orthogonal, 


u m (t) = Z a.u 1 (t) 
i=l 1 


Z a. = 1 
i=l 1 


Using (A. 27), 


-Tr{M* 1 (-2a*H* + Z a.M* 1 )} 

i=l 1 


£ £ 

2a m - 2a m = 0 
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where a is the component of u in u m . Since u is a stationary point of 
log |D| and logarithm is a monotonic function, it must be a stationary 
point for |D| . 


A. 5 LINEAR FUNCTIONAL OF THE DISPERSION MATRIX 

* 

Theorem A. 3 : A necessary condition that an input u (t) minimizes <5?(D) is that 

* * * 

(i) u (t) is an eigenvector of Q with eigenvalue J2?(D ). 

* 

(ii) Any other eigenvalue of Q follows the inequality 


\ 1 (1 + (M* 1 M* k ) 2 )} 


(A. 35) 


and 


Q ± *(t, s) = &(K* V.(t, s)) 


(A. 36) 


Proof : It is similar to proof of Theorem A. 1. 


Theorem A. 4: A sufficient condition that an input u (t) is a stationary point for 

* * * 
j?(D) is that u be an eigenvector of Q with eigenvalue i^(D ) . 


Proof : Again the proof follows along the same lines as in Theorem A. 2, 


A . 6 EXTENSIONS 


We state one more theorem without proof. 


Theorem A. 5: A necessary and sufficient condition than an input u (t) is a 

g * * 

stationary point for J2?(D ) is that u (t) be an eigenvector of with 

*£ * 
eigenvalue j?(D ) , where 
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This is an important theorem because by choosing the linear operator as 
trace and a high value of i , the maximum eigenvalue of D can be minimized. 
In other words , the maximum error in any direction in parameter space is 
minimum. 



APPENDIX B 

INPUT DESIGN IN FREQUENCY DOMAIN 
B . 1 PROBLEM STATEMENT 

Consider the state space representation of a discrete time system 

x(k+l) = $x(k) + Gu(k) k=l,2, s * a N 
and the noisy measurements 

y (k) = Hx(k) + v(k) 

where 

x(*) is n x 1 state vector 
u( * ) is q x 1 control vector 
y ( • ) and v(*) are p x 1 measurement and noise vectors, respectively. 

(p, G, and H are matrices of appropriate dimensions and contain m unknown para- 
meters 0. It is assumed that: 

(a) The system is stable and 0 identifiable. 

(b) v is a second order process with known covariance. 

[ 22 ] 

Mehra developed a technique for determining the spectrum of input u which 
minimizes the trace or determinant of the dispersion matrix of parameter esti- 
mates . A summary is presented here . 


(B . 1) 


(B . 2) 
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B . 2 FREQUENCY DOMAIN APPROACH 


The approach consists of transforming the frequency domain representation 

of the system into a regression problem. The results of Kiefer and Wolfowitz 
[14] 

and Fedorov, in the theory of optimal experiments , are subsequently applied 
to develop schemes for designing optimal inputs . 

Fourier transform (1) and (2) to get: 

e -jtO o n = ^ (n) + G G(n) ( B .3) 


and 


y(n) = Hx(n) + v(n) 

where 

N _ N ,, N , 

n - - 2 ‘, ~2 "r ■*-> • • • » 2 ■ L 

= 2tt/N 

and ~ over a variable denotes its finite Fourier transform. 

From (B.3) and (B.4) 

y(n) = H(e ^ W o n - 4>) * Gu(n) + v(n) 

- T(n,0) u(n) + v(n) (B . 5) 

- J T '(n,0 ) + |~ (0-6 o ) + 0(A9 2 )j u(n) + v(n) 

' • o 

where Q q is the a priori estimate of 0 and u is a scalar. Thus, 
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(B . 6) 


Ay(n) - y(n) - T(n,0 ) u(n) 


= u(n) (0- 0 q ) + v(n) 


From (B .6) we can estimate A0 using generalized least-squares 


A0 = 


1 1 9T* q -l, v'3T ( . 

N Re _ N ^ 2 36 S w (n) 30 S uu (n) 


-1 


N/2-1 * 

rjjj Re l S^Cn) u (°) Ay(n) (B.7) 


and cov(A0) = M where M is the information matrix given by 


N/2-1 * . 

M = Re E S 1 (n) •— S (n) 

-N/2 vv 30 uu 


(B . 8) 


S uu (n) and S vv (n) are the spectra of u and v and * denotes transpose and com- 
plex conjugate. 


The average information matrix per sample is 




N/2-1 

E 

-N/2 


•k 

3T -1 
30 vv 


(n) 


3T 

30 


S (n) 
uu 


(B . 9) 


As the number of sample points increases , the information matrix approaches 
infinity. However, the average information matrix per sample reaches a constant 
value . 


We now take limits of equation (B . 9) as N -*■ 00 . It is straightforward to 
show that Lt M. T -* M , where 

N oo 

N-**> 1 
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(B.10) 





S_1(w) H dF 

w o0 uu 


( 03 ) 


where F is the spectral distribution function of u. In the input design procedure 
the total power in u will be constrained and, without loss of generality, we may 
fix it at unity . Such an input will be called normalized input and the correspond- 
ing information matrix is the normalized information matrix. We will only con- 
sider normalized inputs in the following development. 

The normalized information matrix has many important properties . Some 
of them are given in Theorem 1 (see [2] for proof) . 


Theorem 1: 


(a) M is a symmetric positive semi- definite matrix. 

(b) If the spectrum of u contains fewer than ^ frequencies, the matrix 
M is singular and not all parameters can be identified. 

(c) For any normalized input with mixed spectrum f^(oo) , it is always 
possible to find another normalized input f 2 (oo) which has the same 
information matrix and has discrete spectrum with no more than 

— ^ + 1 frequencies. In addition, if f(oo) is an optimizing nor- 
malized input, the number of frequencies cannot exceed ^ . 


This theorem has important implications . Thus, the optimal input is a sum 
lite number of sinusoidal 
quencies in the optimal design, 


of a finite number of sinusoidal waves . If k is the number of independent fre- 


m , m(nH-l) 
2 P - q- 2 


(B.ll) 


There are two other theorems which are important for input design. 
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Theorem 2; The following are equivalent: 


(a) The normalized design f* maximizes the determinant of the information 
matrix (or minimizes the determinant of the dispersion matrix) . 

(b) The normalized design f* minimizes 

max Tjj(o),f) 
cue [-tt,tt] 

( c ) max ip (to, f *) = m 

U)£[-TT,TT] 

* 

( 1 'Nm _ *1 \ 

5^(03) M ^ ) 

Theorem 3 : The following are equivalent: 


(B . 12) 
(B . 13) 
(B . 14) 


(a) The normalized design f minimizes i2?(D(f*) 

(b) The normalized design f* minimizes 

max ip(oj f f) 

O)0[-1T,TT] 

( c ) ™ax iKu.f*) = i2?(D(f*)) (B 15) 

0)e[-1T,TT] 

where ^(a),f) = (f ) s ^(w) D(f)J (B . 16) 

and 3? is a linear operator (see Equation (3.9) for definition of J?0 . 

All input spectra satisfying (a)-(c) of theorem 2 and their linear combin- 
ations give the same information matrix. The same holds of Theorem 3. 


The above theorems are very important because they convert complex non- 
linear problems into simpler ones. Instead of minimizing |D| or j^TD), it is only 

necessary to minimize the max \J/(co,f) . This leads to a very powerful 

coe [-7t ,71 ] 

computation technique given in Section 3.5. 
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Example 1 : Consider a second order system in which we wish to estimate 

the frequency 0^ 2 



(B . 17) 


y = x i + v 

Discrete time equivalent to this system is (for small A) 



‘ 1 

A" 


"0" 

x(k+l) = 



x(k) + 



_-A0 

-A_ 


_A_ 


u(k) 


(B , 18) 


(B . 19) 


z(k) = x^(k) + v(k) 


(B . 20) 


We assume that v(k) is . white with known power spectral density S 


vv 


T(w,0) = [1 0] 


V*“-l . 0 - 

-1 

’°‘ 

. AS e" 3 “+4. 


_A_ 


(e" jW -l)(e“ jW -l+A) + A 2 9 


(B .21) 



M(w) 


{(e JW -l) e jW -l+A) + A 2 0} 2 
A 8 S" 1 

+ - - XX 

{ (e jW -l) (e jW -l+A)+A 2 0} 2 { (e' jW -l) (e ja) -l+A)+A 2 0} 2 


(B . 22) 


(B.23) 
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In this case, the maximum number of frequencies is one and can be determined 
by maximizing M(co) for coe[-7C,7t] . It can be shown that for small A, M is 
maximized at 





- 0 ^ > 8 > ^ (B .24) 

0 is not smaller than — from the assumption of a stable system. For high damping, 
it is best to use a constant input. With low damping, an oscillatory input at the 
damped natural frequency is the best. 
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APPENDIX C 

COMPUTATION OF TIME DOMAIN INPUTS USING 
EIGENVALUE-EIGENVECTOR DECOMPOSITION 


In the time domain input design problem, a weighted trace of information 
matrix has to be maximized in each iteration, whichever criteria of Section 3.3 
we may be using. This leads to a two point boundary value problem of the type 
(see Equation (3 . 10)) : 

x e 

x • 


r n 


T 1 

x e 


F e 

-yG e G e 

.X _ 


H e E 0 lH e 

1 

CD 

H 

1 


A 



x Q (0) = 0 


X(T) = 0 


(C.l) 


The smallest value of p to give a nontrivial solution to the above equation is to 

[28] 

be determined. It has been shown by Bryson and Hall 1 that the Hamiltonian 
matrix <3$fis symplectic. Therefore, the eigenvalues of occur in pairs +X. and 
-X. Let SP+ and be the positive and negative sets of eigenvalues of Then 

the corresponding eigenvector matrix can be written in the partitioned form as 



L A+ A -J 

The eigenvectors are normalized such that 
A^X + - X^A + = I 


(C . 2) 


(C . 3) 
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It can be shown that (see Bryson arid Hall) 


S 1 = 


a t -x t 


-A 


(C. 4) 


Premultiply both sides of Equation (C . 1) by S 


-1 


T-S- 1 

dt 


~e 

V 


= s"k*ss 1 


(€.5) 


or 


_d_ 

dt 


T T 

A_x 0 - X_A 


T T 

-A + x 0 .+ x;x 


'+ I 


0 I -A 


+J 


a\ q - X^A 


T T 

A + x e + X + X 


CC.6) 


Using the boundary conditions and simplifying, 


C T v e - V T A T 


e- + V '< A I 

+ 1 + 


A(0) 


L x e (T) 


= o 


(C. 7) 


For a nontrivial solution, 


-s , T.T 
e + A 


-s T T 
e + X , 

» T 


is singular 


- e~ S + T X^x‘ T e S + T A* is singular 


i. e. 


, U = A + 1 A_e~ s + T x_ 1 X + e -s + T 


(C . 8) 
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has at least one eigenvalue equal to one. 

Suppose for a certain p and T , U has eigenvalues close to one . Let 

^ ® . f n /'•i\ 

the corresponding normalized left and right eigenvectors of U be J and ' . 

Then 


(1) _ . (1) (1) 


Uy-= X 


(C .9) 


Differentiating with respect to T, 


9 y R 1} 


M y d) + u 

9T y R U 9T 


(1) y^ 1 ) 

U; (1) . ,(1) 9 y R 
9T y R , 9T 


9A 


(C.10) 


Cl)T 

Premultiply by y T and simplify to get 




»«> . /.< D* 

9T \ y L 


Similarly , 


& a) - k i,T 0 ^)k i)T 


(C.ll) 


15 = -A h i.e'Vx’Ve’V - A.h e" i e’ J + T 

dT +-+ — + “T“ * T. T 


9U 3(A + 1 A_) _ s T -! -s.T A -1 A _-s^T 9s + v -l v — s^T,,, 

3iT = ~9TI e +X - X + e + " A + A - e + 3]T X - X + e + * T 

I A -1 A _S I T X+ -T S | T A -1 A n -S | T Y -1 V n -S +^ T 3s+ * T 

+ A + A - e + w 6 + ~ A + - x - x + e + 9ir T 

Since ^ are the first half eigenvalues of 


(C . 12) 
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= A X . -X T . 

\ -i -iy 


0 - G e G e 


iO o 



T T 

- - A -i G e G e A +i 


or 


^ 4 + T T 

air = diag t_A - G 0 G e A + ] 


CC.13) 


It can also be shown that 


where 


9 _ 


X_ 


X + 

X_ 


V 

A- 


A + 

A__ 


(G . 14 ) 


ib . . = - A T . G q Gq A,, /(A. - A.) 1 < j <n 0 , 1 < i<n 0 

ji -J00+1J l — — __ 


■ A la - n0) G e G 0 A +i nx j - V 


= -A T . G-gI A.. a i/(A.-A.) I<j<n 0 , n 0 <i< 2 n 0 

-J 0 9 “(j “ n 0 ) j x - - 


-A T /. qs G„Gq A , . a v /(A. - A.) 

-(j - n 0 ) 00 -(i - n 0 ) j i 


(C . 15) 


n 0 < j <_ 2 n 0 , n 0 < i 2 n 0 
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a(A + 1 A_) acxlX) 

Equation (C , 14) can be used to compute ^ and 




3p 


To find the change in p for a small change AT in T, we use the relation 


9A 


( 1 ) 


9T 


AT + 


9A 


( 1 ) 


3p 


Ay = 0 


(C . 16) 


If the desired eigenvalue is close to one, but not exactly equal to one, the approximate 
change in y. required to bring this eigenvalue closer to one is given by 


9A 


( 1 ) 


9y 


Ay « 1 - 


A (1 > 


(C . 17) 


For a given \i, there may be several values of T for which the matrix U 
has at least one eigenvalue equal to unity. We are interested in the smallest T 
for which this is true. If we start with a correct p^T^ pair, this iteration tech- 
nique will always give the correct p,T pair for small changes, AT. Thus, if the 

desired T is "far" from T , \x is updated in steps, each step involving a small 

° ( 1 ) 

change in T . In each step (C. 17) is used to bring A, closer to one. Once the 
optimum |! is found, A(0) can be determined. Then the state and control time 
histories and information matrix are evaluated. 


This technique gives an excellent insight into the structure of the optimal 
inputs. By studying the eigenvalues and eigenvectors, we can determine if 
they consist of basically damped or undamped oscillatory functions of rising and 
falling exponentials . The time constants of these oscillatory functions and 
exponentials can usually be correlated with model time constants . This approach 
is excellent for obtaining good approximations to optimal inputs and building 
input generators . Instead of specifying the input as a function of time, it is 
possible to give its functional representation. The sensitivity of the nature of 
optimal inputs to the length of the experiment can also be studied. 


Example 


Consider the following first order system with an unknown parameter 0, 
with a priori value one. 
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x = -x + Bu x(0) =0 (C . 18) 

There is a continuous measurement of the scalar state x: 

y = x + v (C . 19) 

Let 

E(v(t)) = 0, E(v(t)v(x) ) = 6(t - t) 




R e = ti] (c . 20 ) 

The controllability matrix of (Fq, Gq) is 

-1 -1 

-1 -1 

The system is uncontrollable. The controllable part can be represented by 
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[- 1 ] 

[ 1 ] 


(C . 21) 


(C . 22) 


normalized eigenvector 


(C.23) 


(C . 24) 


(C .25) 


(C.26) 
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or 


tan 




Equation (C . 2 7 ) has 
falling i n the second 


™any solutions 
quadrant. For 


The smallest T corresponds to 
example, if T = 2, p = 2.30. 


(C .27) 
v /~M~- I T 
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APPENDIX D 


PRACTICAL TECHNIQUES FOR SENSITIVITY FUNCTIONS 
REDUCTION IN LINEAR TIME-INVARIANT SYSTEMS 

D. 1 INTRODUCTION 

The problem of computing state sensitivities using reduced order models 
has become very important in parameter estimation involving high order models 
and many unknown parameters. These techniques allow a considerable saving 

in computation time which makes the determination of optimal inputs feasible for 

[29-311 

practical systems. Most efforts to date have concentrated on finding 

bounds on the order of the model which can generate state sensitivities for all 
system parameters. Very little attention has been given to the formulation of 
practical techniques leading to these lowest order models. Formulations by 
Wilkie and Perkins, Denery, and Neuman and Sood^^ lead to fairly 

complicated transformations and are not capable of exploiting the characteristics 
of the system in most cases . 

A practical method for obtaining lowest order models for sensitivity func- 
tions computations is developed. The technique makes full use of special system 
characteristics and has general application to high order systems with a large 
number of unknown parameters . 

D . 2 PROBLEM STATEMENT 

Consider a system 

x = Fx + Gu x(0) - x (D.l) 

o 

where x is a nxl state vector, u is a qxl control vector. F and G are matrices 
of appropriate dimensions and are functions of m parameters 0. 
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A heretofore uncited property of systems , which depends on the para- 
meters 0, is important in sensitivity computation. 

Definition 1 - Structural Controllability : A system is said to be structurally 

controllable if it is controllable for almost all values of parameters . The 
system may be uncontrollable if certain relations hold among the para- 
meters . 


Definition 2 - Structural Linear Dependence : A set of vectors has structural 
linear dependence if a linear combination of these vectors is zero for 
almost all values of parameters . The particular linear combination may 
depend on the values of the parameters . 


Ex amt 


Consider the system 



(D . 2) 


The controllability matrix is 



The system is controllable unless 0^ + 0^ = 0^ + 0^. Thus, if 0^ = 0^ = -1 and 
0^ = 0^ = -5, . the system is uncontrollable in the classical sense but structurally 
controllable. 

Initially, the following simplifications can be made: 

a. The system is made structurally controllable by dropping uncon- 
trollable states. Since the initial condition is zero, the system 
never moves into the uncontrollable subspace, This reduces the 
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order of the system. Note that the states which are uncontrollable 
only for the given values of the parameters but which are structurally 
controllable should not be dropped. 

b . All structurally linearly dependent columns of G matrix are lumped 
with other columns. This reduces the number of effective controls. 

[29] 

The state sensitivities for all parameters 0 can be written as 


x e = F e x e + V 


(D . 3) 


x 0 (O) = 0 

Fg , Gg, and Xg are defined in Chapter III. If (Fq, Gg) is uncontrollable, the 
corresponding controllability matrix is of rank less than (m+l)n, say r. Let 
Q-^ be the set of r independent columns in the controllability matrix. If 
is such that and Q ^ form a set of n(m + l) linearly independent vectors, then 


/-t\ 


x 6 4 <«! i a 


V Q 2/ 


follows the differential equation 


, 

F n 

\° 



(D . 4) 


(D . 5) 


Since the initial condition in (D.5) is zero, the last (m+l)n-r uncontrollable 
states remain zero throughout. The remaining r states, x^, follow the differential 
equation 
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(D .6) 


X c = F o * + G u 

c c c c 


X (0) = o 

c 


where 

F c - F n = Q i F e Q i 

(D . 7) 

g c a G ; = q }g 0 

Also, 


x 8 ■ «i ; v * 8 


Ve 


(D.8) 


since other states in x 0 are zero. Note that is a pseudo-inverse of 
depending on Q 2 . The transformation from F 0 ,G 0 to F c ,G c and from x £ to 
Xq does not involve Qg explicitly . Therefore , can be chosen to be any 
pseudo-inverse of Q^, for example, 


Q 



(Q* Q 1 )" 1 Q* 


(D. 9) 


It is assumed here that the inputs are linearly independent. If this is not 
so, the number of inputs can be reduced until they are linearly independent. 

This will usually result in a reduction in the controllable subspace of (F 0 ,G 0 ) 
as shown in Section D .4. Under the assumption of linear independence of inputs, 
it is necessary and sufficient to solve a system of r linear equations (D.6) to 
determine the state and its sensitivities at all times. The next sections investi- 
gate the nature and dimension of the controllable subspace of (F 0 ,G 0 ) and 
explore efficient methods for finding Q^. 
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D . 3 SINGLE INPUT SYSTEMS 


In single input systems, Fq is a (m+ln x (m+l)n matrix and Gq is a 
(m+l)n vector . The following theorem holds . 

Theorem D . 1 : For a single input system, the rank of the controllability matrix 

of (FqjCjq) is less than 2n. 

Proof : The controllability matrix of (Fq,Gq) is 


C. - 


G e : F o c 0 ; 


It is easy to show that 


F P G 


If 


F e G e ' 


2-f p g 


90 


9 „p 

3T FG 

m 


F (m + l)n - 1 

0 G 0 


A D(F P G) A + D * (p P G) 


(D.10) 


(D.ll) 


T^n 

F = a I + a. F < 
o 1 


a ,F 
n-1 


n-1 


(D . 12) 


the (n+k)th column of Cq is 


m+k-1 


n-1 


D(F 11 ^ J 'G) = Z {D"(a DF^g + a. D (F lC+1 ’"’ 1 G) } 
i=0 1 1 


ck+i-1 


(D . 13) 
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Thus , the order of the system required to compute all state sensitivities 
for a single input system cannot exceed 2n. In many practical cases, it is smaller 
as shown in Example 1. 

Corollary 1 : If the structurally controllable subspace of (F,G) is of the order 

p, the maximal order of the controllable subspace of (FqjGq) is 2p. 

Since F^G is a linear combination of G, FG, . . . , F^ ^G, the corollary follows 
immediately . 
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Example 2 : Consider the following system 


• / 9 1 0 2 ' 
x =( ^ Z |x +| 

\0 -] 


U 


x(0) = 0 


(D . 17) 


The state vector and its sensitivities for 0-^ and 0£ form a set of six differential 
equations . Since the number of states is two and the number of controls is one , 
only the first four columns can be independent in the controllability matrix of 
(F 0> G 0 ) . These columns are 


rank(C ) = rank 
8 


»(;). °Q). ' 82 ). "( 


e i 0 2 - 6 i e 2 

-1 


-2) 


(D . 18) 


rank 


1 

0 

0 

0 


-1 

0 

0 

1 


"A - e 2 


1 

0 , 


e i-i 


e i 0 2 - e ! e 2 + 6 2 


-1 


20 0 - 6 
1 2 2 


9 i - e i + 1 


(D. 19) 


The first three columns are independent for ©2 ^ 0* If ©2 zero > only the first 
two columns are independent and the required model is of the second order . 


D . 4 MULTI-INPUT SYSTEMS 


We state and prove the following theorem . 

Theorem D , 2 ; The rank of the controllability matrix of (Fq,Gq) cannot exceed 
(q+l)n for a q inputs system. 
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In multi-input systems, Gq is a matrix with q columns. The controllability 
matrix of (F e* G e ) has (m+l)nq columns which can be written thus: 


= [D(G 1 ,FG 1> . 


>F n<m+l) 1 Gi ):d(C 2 ,FC 2> . 


pnO'rt-D-lgj,. 


...:D(G ,FG 
q q 


<tF n(m+l)-l }] 


(D . 20) 


where is the ith column of G. By using Theorem 1, it is easily seen that the 
last (m-l)n columns involving any of vectors G/s are linearly dependent on the 
first 2n columns for that vector. From (D.14) and (D . 15) 


rank(Cg) - p = rank[{D(G 1 ,FG 1 F n ~ 1 G 1 ),ZD*(a i I)F 1 G 1 ZD* (a i I)F 1+n 1 G 1 > 

,..{D(G ,FG ,...,F n-1 G ),ZD*(a.I)F i G n> . .. ,ZD*(a.I)F 1+n ~ 1 G» 

q q q i q 1 q 

(d . 21) 


where all summations are from 0 to n-1. Let the structuredly independent columns 
in the controllability matrix of (F,G) be 


[G 


1’ 


,F k l" 1 G 


r 


.,F k 2 _1 G 


2* 




* F k q -1 G q ] 


(D . 22) 


q 

Z k. = n 
i=l 1 


(D . 23) 


This set of linearly independent vectors in the controllability matrix spans the 
complete n- dimensional space . So any vector can be represented as a linear 
combination of these vectors . In particular, 



3 1 G 1 + ’ 


.+ s k * k i 

k l 



6 F k q X G 
n q 


1 < j < n 
1 < k < q 


(D . 24) 
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Therefore, 


n-1 .... n-1 . n-1 

2 D*(a I)F 1+J V, = 6. E D* (a. I)F 1 G. + B, Z D*(a.I)F 1+1 G, 

1=0 1 i =0 1 1 Z i =0 1 1 


n-1 

, + B E D*(a.I)F q G 
n 1=0 1 « 


(D . 25) 


This is a linear combination of n vectors in the right hand side matrix of (D.24) 
for all j and k (the values of p^ depend on j and k) . Thus, 


p rank[D(G^,FG^, . 


,F n 1 g 1 ):d(g 2 ,fg 2 , 


. . ,F n-1 G„) : . . . . :d(g ,FG 

2 . • q q 



ZD* (a I) F 1 G 1 , ED* (c^i) f 1+1 g 


,ED* (a. I)F G. 
i 3 


< (q+l)n 


,ED*(a i I)F 1+k q X G ] 
(D .26) 


Thus , the procedure for finding independent columns of Cg consists of 
finding structurally independent columns of the controllability matrix of (F,G) 
choosing (q+l)n appropriate columns from Cg and checking to see if there is any 
further linear dependence. 


Another simplification is possible in large order systems in which each 
input controls only a small number of states . If k! is the dimension of the con- 
trollable subspace of the ith input, no more than 2ki columns involving can 
be linearly independent in the right hand side matrix of (D.26) as shown in 
Corollary 1 to Theorem 1. 

Corollary 2 : If for any single input u^ the system is completely controllable , 

p = rank[D(G 1# . . . ,F n-1 G ;L ) ,D(G 2> . . . ,F n-1 G 2 ) . . ,D(G , . . . ,F n_1 G q ) , 

.... ,ED*(a.I)F i G j , . . . ,ED*(a.I)F nH ' 1 ' 1 G j ] (D .27) 
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The proof is obvious since G., . . 
pendent columns . 



Gj form a set of n linearly inde- 


D . 5 CONCLUSIONS 


A systematic method for finding the controllable subspace of the augmented 
system, in which the state vector is the system state and its sensitivities, is 
presented. It is necessary to start with no more than (q+l)n columns of the 
controllability matrix of the augmented system . These columns can be selected 
quickly by inspection of the controllability matrix of the initial system. If r 
is the dimension of the controllable subspace of the augmented system, it is 
necessary to solve r linear equations to evaluate the state vector and its sensitivities . 

This method of sensitivity functions reduction fully exploits the character- 
istics of the system and the cases in which the sensitivity to all parameters in 
the system is not required. In other words, it leads to the minimal order model 
under the circumstances. 
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